/[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 121 by jgs, Fri May 6 04:26:16 2005 UTC trunk/escript/src/Data.cpp revision 1872 by jfenwick, Mon Oct 13 00:18:55 2008 UTC
# Line 1  Line 1 
 // $Id$  
 /*=============================================================================  
1    
2   ******************************************************************************  /*******************************************************
3   *                                                                            *  *
4   *       COPYRIGHT ACcESS 2004 -  All Rights Reserved                         *  * Copyright (c) 2003-2008 by University of Queensland
5   *                                                                            *  * Earth Systems Science Computational Center (ESSCC)
6   * This software is the property of ACcESS.  No part of this code             *  * http://www.uq.edu.au/esscc
7   * may be copied in any form or by any means without the expressed written    *  *
8   * consent of ACcESS.  Copying, use or modification of this software          *  * Primary Business: Queensland, Australia
9   * by any unauthorised person is illegal unless that                          *  * Licensed under the Open Software License version 3.0
10   * person has a software license agreement with ACcESS.                       *  * http://www.opensource.org/licenses/osl-3.0.php
11   *                                                                            *  *
12   ******************************************************************************  *******************************************************/
13    
14    
15    #include "Data.h"
16    
17    #include "DataExpanded.h"
18    #include "DataConstant.h"
19    #include "DataTagged.h"
20    #include "DataEmpty.h"
21    #include "FunctionSpaceFactory.h"
22    #include "AbstractContinuousDomain.h"
23    #include "UnaryFuncs.h"
24    #include "FunctionSpaceException.h"
25    
26  ******************************************************************************/  extern "C" {
27    #include "escript/blocktimer.h"
28  #include "escript/Data/Data.h"  }
29    
 #include <iostream>  
30  #include <fstream>  #include <fstream>
31  #include <algorithm>  #include <algorithm>
32  #include <vector>  #include <vector>
 #include <exception>  
33  #include <functional>  #include <functional>
 #include <math.h>  
34    
35  #include <boost/python/str.hpp>  #include <boost/python/dict.hpp>
36  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
37  #include <boost/python/long.hpp>  #include <boost/python/long.hpp>
38    
 #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"  
   
39  using namespace std;  using namespace std;
40  using namespace boost::python;  using namespace boost::python;
41  using namespace boost;  using namespace boost;
# Line 50  Data::Data() Line 46  Data::Data()
46    //    //
47    // Default data is type DataEmpty    // Default data is type DataEmpty
48    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
49    shared_ptr<DataAbstract> temp_data(temp);    m_data=temp->getPtr();
50    m_data=temp_data;    m_protected=false;
51  }  }
52    
53  Data::Data(double value,  Data::Data(double value,
# Line 59  Data::Data(double value, Line 55  Data::Data(double value,
55             const FunctionSpace& what,             const FunctionSpace& what,
56             bool expanded)             bool expanded)
57  {  {
58    DataArrayView::ShapeType dataPointShape;    DataTypes::ShapeType dataPointShape;
59    for (int i = 0; i < shape.attr("__len__")(); ++i) {    for (int i = 0; i < shape.attr("__len__")(); ++i) {
60      dataPointShape.push_back(extract<const int>(shape[i]));      dataPointShape.push_back(extract<const int>(shape[i]));
61    }    }
62    DataArray temp(dataPointShape,value);  
63    initialise(temp.getView(),what,expanded);    int len = DataTypes::noValues(dataPointShape);
64      DataVector temp_data(len,value,len);
65      initialise(temp_data, dataPointShape, what, expanded);
66      m_protected=false;
67  }  }
68    
69  Data::Data(double value,  Data::Data(double value,
70         const DataArrayView::ShapeType& dataPointShape,         const DataTypes::ShapeType& dataPointShape,
71         const FunctionSpace& what,         const FunctionSpace& what,
72             bool expanded)             bool expanded)
73  {  {
74    DataArray temp(dataPointShape,value);    int len = DataTypes::noValues(dataPointShape);
75    pair<int,int> dataShape=what.getDataShape();  
76    initialise(temp.getView(),what,expanded);    DataVector temp_data(len,value,len);
77    //   DataArrayView temp_dataView(temp_data, dataPointShape);
78    
79    //   initialise(temp_dataView, what, expanded);
80      initialise(temp_data, dataPointShape, what, expanded);
81    
82      m_protected=false;
83  }  }
84    
85  Data::Data(const Data& inData)  Data::Data(const Data& inData)
86  {  {
87    m_data=inData.m_data;    m_data=inData.m_data;
88      m_protected=inData.isProtected();
89  }  }
90    
91    
92  Data::Data(const Data& inData,  Data::Data(const Data& inData,
93             const DataArrayView::RegionType& region)             const DataTypes::RegionType& region)
94  {  {
95    //    //
96    // Create Data which is a slice of another Data    // Create Data which is a slice of another Data
97    DataAbstract* tmp = inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
98    shared_ptr<DataAbstract> temp_data(tmp);    m_data=DataAbstract_ptr(tmp);
99    m_data=temp_data;    m_protected=false;
100  }  }
101    
102  Data::Data(const Data& inData,  Data::Data(const Data& inData,
103             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
104  {  {
105      if (inData.isEmpty())
106      {
107        throw DataException("Error - will not interpolate for instances of DataEmpty.");
108      }
109    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
110      m_data=inData.m_data;      m_data=inData.m_data;
111      } else if (inData.isConstant()) { // for a constant function, we just need to use the new function space
112        if (!inData.probeInterpolation(functionspace))
113        {           // Even though this is constant, we still need to check whether interpolation is allowed
114        throw FunctionSpaceException("Call to probeInterpolation returned false for DataConstant.");
115        }
116        DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),inData.m_data->getVector());  
117        m_data=DataAbstract_ptr(dc);
118    } else {    } else {
119      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getDataPointShape(),functionspace,true);
120      // Note for Lutz, Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
121      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
122      // be able to create an instance of AbstractDomain but that was done      // be able to create an instance of AbstractDomain but that was done
123      // as a boost python work around which may no longer be required.      // as a boost:python work around which may no longer be required.
124      const AbstractDomain& inDataDomain=inData.getDomain();      /*const AbstractDomain& inDataDomain=inData.getDomain();*/
125        const_Domain_ptr inDataDomain=inData.getDomain();
126      if  (inDataDomain==functionspace.getDomain()) {      if  (inDataDomain==functionspace.getDomain()) {
127        inDataDomain.interpolateOnDomain(tmp,inData);        inDataDomain->interpolateOnDomain(tmp,inData);
128      } else {      } else {
129        inDataDomain.interpolateACross(tmp,inData);        inDataDomain->interpolateACross(tmp,inData);
130      }      }
131      m_data=tmp.m_data;      m_data=tmp.m_data;
132    }    }
133      m_protected=false;
134  }  }
135    
136  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(DataAbstract* underlyingdata)
            const DataTagged::ValueListType & values,  
            const DataArrayView& defaultValue,  
            const FunctionSpace& what,  
            bool expanded)  
137  {  {
138    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);  //  m_data=shared_ptr<DataAbstract>(underlyingdata);
139    shared_ptr<DataAbstract> temp_data(temp);      m_data=underlyingdata->getPtr();
140    m_data=temp_data;      m_protected=false;
   if (expanded) {  
     expand();  
   }  
141  }  }
142    
143  Data::Data(const numeric::array& value,  Data::Data(const numeric::array& value,
# Line 132  Data::Data(const numeric::array& value, Line 145  Data::Data(const numeric::array& value,
145             bool expanded)             bool expanded)
146  {  {
147    initialise(value,what,expanded);    initialise(value,what,expanded);
148      m_protected=false;
149  }  }
150    /*
151  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
152         const FunctionSpace& what,         const FunctionSpace& what,
153             bool expanded)             bool expanded)
154  {  {
155    initialise(value,what,expanded);    initialise(value,what,expanded);
156      m_protected=false;
157    }*/
158    
159    Data::Data(const DataTypes::ValueType& value,
160             const DataTypes::ShapeType& shape,
161                     const FunctionSpace& what,
162                     bool expanded)
163    {
164       initialise(value,shape,what,expanded);
165       m_protected=false;
166  }  }
167    
168    
169  Data::Data(const object& value,  Data::Data(const object& value,
170         const FunctionSpace& what,         const FunctionSpace& what,
171             bool expanded)             bool expanded)
172  {  {
173    numeric::array asNumArray(value);    numeric::array asNumArray(value);
174    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
175      m_protected=false;
176  }  }
177    
178    
179  Data::Data(const object& value,  Data::Data(const object& value,
180             const Data& other)             const Data& other)
181  {  {
182      numeric::array asNumArray(value);
183    
184      // extract the shape of the numarray
185      DataTypes::ShapeType tempShape=DataTypes::shapeFromNumArray(asNumArray);
186    // /*  for (int i=0; i < asNumArray.getrank(); i++) {
187    //     tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
188    //   }*/
189    //   // get the space for the data vector
190    //   int len = DataTypes::noValues(tempShape);
191    //   DataVector temp_data(len, 0.0, len);
192    // /*  DataArrayView temp_dataView(temp_data, tempShape);
193    //   temp_dataView.copy(asNumArray);*/
194    //   temp_data.copyFromNumArray(asNumArray);
195    
196    //    //
197    // Create DataConstant using the given value and all other parameters    // Create DataConstant using the given value and all other parameters
198    // copied from other. If value is a rank 0 object this Data    // copied from other. If value is a rank 0 object this Data
199    // will assume the point data shape of other.    // will assume the point data shape of other.
200    DataArray temp(value);  
201    if (temp.getView().getRank()==0) {    if (DataTypes::getRank(tempShape)/*temp_dataView.getRank()*/==0) {
202      //  
203      // Create a DataArray with the scalar value for all elements  
204      DataArray temp2(other.getPointDataView().getShape(),temp.getView()());      // get the space for the data vector
205      initialise(temp2.getView(),other.getFunctionSpace(),false);      int len1 = DataTypes::noValues(tempShape);
206        DataVector temp_data(len1, 0.0, len1);
207        temp_data.copyFromNumArray(asNumArray);
208    
209        int len = DataTypes::noValues(other.getDataPointShape());
210    
211        DataVector temp2_data(len, temp_data[0]/*temp_dataView()*/, len);
212        //DataArrayView temp2_dataView(temp2_data, other.getPointDataView().getShape());
213    //     initialise(temp2_dataView, other.getFunctionSpace(), false);
214    
215        DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
216    //     boost::shared_ptr<DataAbstract> sp(t);
217    //     m_data=sp;
218        m_data=DataAbstract_ptr(t);
219    
220    } else {    } else {
221      //      //
222      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
223      initialise(temp.getView(),other.getFunctionSpace(),false);  //     initialise(temp_dataView, other.getFunctionSpace(), false);
224        DataConstant* t=new DataConstant(asNumArray,other.getFunctionSpace());
225    //     boost::shared_ptr<DataAbstract> sp(t);
226    //     m_data=sp;
227        m_data=DataAbstract_ptr(t);
228    }    }
229      m_protected=false;
230  }  }
231    
232    Data::~Data()
233    {
234    
235    }
236    
237    
238    
239    void
240    Data::initialise(const boost::python::numeric::array& value,
241                     const FunctionSpace& what,
242                     bool expanded)
243    {
244      //
245      // Construct a Data object of the appropriate type.
246      // Construct the object first as there seems to be a bug which causes
247      // undefined behaviour if an exception is thrown during construction
248      // within the shared_ptr constructor.
249      if (expanded) {
250        DataAbstract* temp=new DataExpanded(value, what);
251    //     boost::shared_ptr<DataAbstract> temp_data(temp);
252    //     m_data=temp_data;
253        m_data=temp->getPtr();
254      } else {
255        DataAbstract* temp=new DataConstant(value, what);
256    //     boost::shared_ptr<DataAbstract> temp_data(temp);
257    //     m_data=temp_data;
258        m_data=temp->getPtr();
259      }
260    }
261    
262    
263    void
264    Data::initialise(const DataTypes::ValueType& value,
265             const DataTypes::ShapeType& shape,
266                     const FunctionSpace& what,
267                     bool expanded)
268    {
269      //
270      // Construct a Data object of the appropriate type.
271      // Construct the object first as there seems to be a bug which causes
272      // undefined behaviour if an exception is thrown during construction
273      // within the shared_ptr constructor.
274      if (expanded) {
275        DataAbstract* temp=new DataExpanded(what, shape, value);
276    //     boost::shared_ptr<DataAbstract> temp_data(temp);
277    //     m_data=temp_data;
278        m_data=temp->getPtr();
279      } else {
280        DataAbstract* temp=new DataConstant(what, shape, value);
281    //     boost::shared_ptr<DataAbstract> temp_data(temp);
282    //     m_data=temp_data;
283        m_data=temp->getPtr();
284      }
285    }
286    
287    
288    // void
289    // Data::CompareDebug(const Data& rd)
290    // {
291    //  using namespace std;
292    //  bool mismatch=false;
293    //  std::cout << "Comparing left and right" << endl;
294    //  const DataTagged* left=dynamic_cast<DataTagged*>(m_data.get());
295    //  const DataTagged* right=dynamic_cast<DataTagged*>(rd.m_data.get());
296    //  
297    //  if (left==0)
298    //  {
299    //      cout << "left arg is not a DataTagged\n";
300    //      return;
301    //  }
302    //  
303    //  if (right==0)
304    //  {
305    //      cout << "right arg is not a DataTagged\n";
306    //      return;
307    //  }
308    //  cout << "Num elements=" << left->getVector().size() << ":" << right->getVector().size() << std::endl;
309    //  cout << "Shapes ";
310    //  if (left->getShape()==right->getShape())
311    //  {
312    //      cout << "ok\n";
313    //  }
314    //  else
315    //  {
316    //      cout << "Problem: shapes do not match\n";
317    //      mismatch=true;
318    //  }
319    //  int lim=left->getVector().size();
320    //  if (right->getVector().size()) lim=right->getVector().size();
321    //  for (int i=0;i<lim;++i)
322    //  {
323    //      if (left->getVector()[i]!=right->getVector()[i])
324    //      {
325    //          cout << "[" << i << "] value mismatch " << left->getVector()[i] << ":" << right->getVector()[i] << endl;
326    //          mismatch=true;
327    //      }
328    //  }
329    //
330    //  // still need to check the tag map
331    //  // also need to watch what is happening to function spaces, are they copied or what?
332    //
333    //  const DataTagged::DataMapType& mapleft=left->getTagLookup();
334    //  const DataTagged::DataMapType& mapright=right->getTagLookup();
335    //
336    //  if (mapleft.size()!=mapright.size())
337    //  {
338    //      cout << "Maps are different sizes " << mapleft.size() << ":" << mapright.size() << endl;
339    //      mismatch=true;
340    //      cout << "Left map\n";
341    //      DataTagged::DataMapType::const_iterator i,j;
342    //      for (i=mapleft.begin();i!=mapleft.end();++i) {
343    //          cout << "(" << i->first << "=>" << i->second << ")\n";
344    //      }
345    //      cout << "Right map\n";
346    //      for (i=mapright.begin();i!=mapright.end();++i) {
347    //          cout << "(" << i->first << "=>" << i->second << ")\n";
348    //      }
349    //      cout << "End map\n";
350    //
351    //  }
352    //
353    //  DataTagged::DataMapType::const_iterator i,j;
354    //  for (i=mapleft.begin(),j=mapright.begin();i!=mapleft.end() && j!=mapright.end();++i,++j) {
355    //     if ((i->first!=j->first) || (i->second!=j->second))
356    //     {
357    //      cout << "(" << i->first << "=>" << i->second << ")";
358    //      cout << ":(" << j->first << "=>" << j->second << ") ";
359    //      mismatch=true;
360    //            }
361    //  }
362    //  if (mismatch)
363    //  {
364    //      cout << "#Mismatch\n";
365    //  }
366    // }
367    
368  escriptDataC  escriptDataC
369  Data::getDataC()  Data::getDataC()
370  {  {
# Line 188  Data::getDataC() const Line 384  Data::getDataC() const
384  const boost::python::tuple  const boost::python::tuple
385  Data::getShapeTuple() const  Data::getShapeTuple() const
386  {  {
387    const DataArrayView::ShapeType& shape=getDataPointShape();    const DataTypes::ShapeType& shape=getDataPointShape();
388    switch(getDataPointRank()) {    switch(getDataPointRank()) {
389       case 0:       case 0:
390          return make_tuple();          return make_tuple();
# Line 205  Data::getShapeTuple() const Line 401  Data::getShapeTuple() const
401    }    }
402  }  }
403    
404    
405    // The different name is needed because boost has trouble with overloaded functions.
406    // It can't work out what type the function is based soley on its name.
407    // There are ways to fix this involving creating function pointer variables for each form
408    // but there doesn't seem to be a need given that the methods have the same name from the python point of view
409    Data*
410    Data::copySelf()
411    {
412       DataAbstract* temp=m_data->deepCopy();
413       return new Data(temp);
414    }
415    
416  void  void
417  Data::copy(const Data& other)  Data::copy(const Data& other)
418  {  {
419    //    DataAbstract* temp=other.m_data->deepCopy();
420    // Perform a deep copy    DataAbstract_ptr p=temp->getPtr();
421      m_data=p;
422    }
423    
424    
425    void
426    Data::setToZero()
427    {
428      if (isEmpty())
429    {    {
430      DataExpanded* temp=dynamic_cast<DataExpanded*>(other.m_data.get());       throw DataException("Error - Operations not permitted on instances of DataEmpty.");
     if (temp!=0) {  
       //  
       // Construct a DataExpanded copy  
       DataAbstract* newData=new DataExpanded(*temp);  
       shared_ptr<DataAbstract> temp_data(newData);  
       m_data=temp_data;  
       return;  
     }  
431    }    }
432    {    {
433      DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());      DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
434      if (temp!=0) {      if (temp!=0) {
435        //         temp->setToZero();
436        // Construct a DataTagged copy         return;
       DataAbstract* newData=new DataTagged(*temp);  
       shared_ptr<DataAbstract> temp_data(newData);  
       m_data=temp_data;  
       return;  
437      }      }
438    }    }
439    {    {
440      DataConstant* temp=dynamic_cast<DataConstant*>(other.m_data.get());      DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
441      if (temp!=0) {      if (temp!=0) {
442        //        temp->setToZero();
       // Construct a DataConstant copy  
       DataAbstract* newData=new DataConstant(*temp);  
       shared_ptr<DataAbstract> temp_data(newData);  
       m_data=temp_data;  
443        return;        return;
444      }      }
445    }    }
446    {    {
447      DataEmpty* temp=dynamic_cast<DataEmpty*>(other.m_data.get());      DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
448      if (temp!=0) {      if (temp!=0) {
449        //        temp->setToZero();
       // Construct a DataEmpty copy  
       DataAbstract* newData=new DataEmpty();  
       shared_ptr<DataAbstract> temp_data(newData);  
       m_data=temp_data;  
450        return;        return;
451      }      }
452    }    }
453    throw DataException("Error - Copy not implemented for this Data type.");    throw DataException("Error - Data can not be set to zero.");
454  }  }
455    
456    // void
457    // Data::copyWithMask(const Data& other,
458    //                    const Data& mask)
459    // {
460    //   if (other.isEmpty() || mask.isEmpty())
461    //   {
462    //  throw DataException("Error - copyWithMask not permitted using instances of DataEmpty.");
463    //   }
464    //   Data mask1;
465    //   Data mask2;
466    //   mask1 = mask.wherePositive();
467    //
468    //   mask2.copy(mask1);
469    //   mask1 *= other;
470    //
471    //   mask2 *= *this;
472    //   mask2 = *this - mask2;
473    //   *this = mask1 + mask2;
474    // }
475    
476  void  void
477  Data::copyWithMask(const Data& other,  Data::copyWithMask(const Data& other,
478                     const Data& mask)                     const Data& mask)
479  {  {
480    Data mask1;    // 1. Interpolate if required so all Datas use the same FS as this
481    Data mask2;    // 2. Tag or Expand so that all Data's are the same type
482      // 3. Iterate over the data vectors copying values where mask is >0
483    mask1 = mask.wherePositive();    if (other.isEmpty() || mask.isEmpty())
484    mask2.copy(mask1);    {
485        throw DataException("Error - copyWithMask not permitted using instances of DataEmpty.");
486      }
487      Data other2(other);
488      Data mask2(mask);
489      FunctionSpace myFS=getFunctionSpace();
490      FunctionSpace oFS=other2.getFunctionSpace();
491      FunctionSpace mFS=mask2.getFunctionSpace();
492      if (oFS!=myFS)
493      {
494         if (other2.probeInterpolation(myFS))
495         {
496        other2=other2.interpolate(myFS);
497         }
498         else
499         {
500        throw DataException("Error - copyWithMask: other FunctionSpace is not compatible with this one.");
501         }
502      }
503      if (mFS!=myFS)
504      {
505         if (mask2.probeInterpolation(myFS))
506         {
507        mask2=mask2.interpolate(myFS);
508         }
509         else
510         {
511        throw DataException("Error - copyWithMask: mask FunctionSpace is not compatible with this one.");
512         }
513      }
514                // Ensure that all args have the same type
515      if (this->isExpanded() || mask2.isExpanded() || other2.isExpanded())
516      {
517        this->expand();
518        other2.expand();
519        mask2.expand();
520      }
521      else if (this->isTagged() || mask2.isTagged() || other2.isTagged())
522      {
523        this->tag();
524        other2.tag();
525        mask2.tag();
526      }
527      else if (this->isConstant() && mask2.isConstant() && other2.isConstant())
528      {
529      }
530      else
531      {
532        throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
533      }
534      // Now we iterate over the elements
535      DataVector& self=m_data->getVector();
536      const DataVector& ovec=other2.m_data->getVector();
537      const DataVector& mvec=mask2.m_data->getVector();
538      if ((self.size()!=ovec.size()) || (self.size()!=mvec.size()))
539      {
540        throw DataException("Error - size mismatch in arguments to copyWithMask.");
541      }
542      size_t num_points=self.size();
543      long i;
544      #pragma omp parallel for private(i) schedule(static)
545      for (i=0;i<num_points;++i)
546      {
547        if (mvec[i]>0)
548        {
549           self[i]=ovec[i];
550        }
551      }
552    }
553    
   mask1 *= other;  
   mask2 *= *this;  
   mask2 = *this - mask2;  
554    
   *this = mask1 + mask2;  
 }  
555    
556  bool  bool
557  Data::isExpanded() const  Data::isExpanded() const
# Line 303  Data::isConstant() const Line 582  Data::isConstant() const
582  }  }
583    
584  void  void
585    Data::setProtection()
586    {
587       m_protected=true;
588    }
589    
590    bool
591    Data::isProtected() const
592    {
593       return m_protected;
594    }
595    
596    
597    
598    void
599  Data::expand()  Data::expand()
600  {  {
601    if (isConstant()) {    if (isConstant()) {
602      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
603      DataAbstract* temp=new DataExpanded(*tempDataConst);      DataAbstract* temp=new DataExpanded(*tempDataConst);
604      shared_ptr<DataAbstract> temp_data(temp);  //     shared_ptr<DataAbstract> temp_data(temp);
605      m_data=temp_data;  //     m_data=temp_data;
606        m_data=temp->getPtr();
607    } else if (isTagged()) {    } else if (isTagged()) {
608      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
609      DataAbstract* temp=new DataExpanded(*tempDataTag);      DataAbstract* temp=new DataExpanded(*tempDataTag);
610      shared_ptr<DataAbstract> temp_data(temp);  //     shared_ptr<DataAbstract> temp_data(temp);
611      m_data=temp_data;  //     m_data=temp_data;
612        m_data=temp->getPtr();
613    } else if (isExpanded()) {    } else if (isExpanded()) {
614      //      //
615      // do nothing      // do nothing
# Line 331  Data::tag() Line 626  Data::tag()
626    if (isConstant()) {    if (isConstant()) {
627      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
628      DataAbstract* temp=new DataTagged(*tempDataConst);      DataAbstract* temp=new DataTagged(*tempDataConst);
629      shared_ptr<DataAbstract> temp_data(temp);  //     shared_ptr<DataAbstract> temp_data(temp);
630      m_data=temp_data;  //     m_data=temp_data;
631        m_data=temp->getPtr();
632    } else if (isTagged()) {    } else if (isTagged()) {
633      // do nothing      // do nothing
634    } else if (isExpanded()) {    } else if (isExpanded()) {
# Line 344  Data::tag() Line 640  Data::tag()
640    }    }
641  }  }
642    
643  void  Data
644  Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  Data::oneOver() const
645  {  {
646    m_data->reshapeDataPoint(shape);    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
647  }  }
648    
649  Data  Data
650  Data::wherePositive() const  Data::wherePositive() const
651  {  {
652    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
653  }  }
654    
655  Data  Data
656  Data::whereNegative() const  Data::whereNegative() const
657  {  {
658    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
659  }  }
660    
661  Data  Data
662  Data::whereNonNegative() const  Data::whereNonNegative() const
663  {  {
664    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
665  }  }
666    
667  Data  Data
668  Data::whereNonPositive() const  Data::whereNonPositive() const
669  {  {
670    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
671  }  }
672    
673  Data  Data
674  Data::whereZero() const  Data::whereZero(double tol) const
675  {  {
676    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));    Data dataAbs=abs();
677      return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
678  }  }
679    
680  Data  Data
681  Data::whereNonZero() const  Data::whereNonZero(double tol) const
682  {  {
683    return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));    Data dataAbs=abs();
684      return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
685  }  }
686    
687  Data  Data
# Line 398  Data::probeInterpolation(const FunctionS Line 696  Data::probeInterpolation(const FunctionS
696    if (getFunctionSpace()==functionspace) {    if (getFunctionSpace()==functionspace) {
697      return true;      return true;
698    } else {    } else {
699      const AbstractDomain& domain=getDomain();      const_Domain_ptr domain=getDomain();
700      if  (domain==functionspace.getDomain()) {      if  (*domain==*functionspace.getDomain()) {
701        return domain.probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());        return domain->probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());
702      } else {      } else {
703        return domain.probeInterpolationACross(getFunctionSpace().getTypeCode(),functionspace.getDomain(),functionspace.getTypeCode());        return domain->probeInterpolationACross(getFunctionSpace().getTypeCode(),*(functionspace.getDomain()),functionspace.getTypeCode());
704      }      }
705    }    }
706  }  }
# Line 410  Data::probeInterpolation(const FunctionS Line 708  Data::probeInterpolation(const FunctionS
708  Data  Data
709  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
710  {  {
711      if (isEmpty())
712      {
713        throw DataException("Error - operation not permitted on instances of DataEmpty.");
714      }
715      double blocktimer_start = blocktimer_time();
716    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
717      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
718    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();    DataTypes::ShapeType grad_shape=getDataPointShape();
719    grad_shape.push_back(functionspace.getDim());    grad_shape.push_back(functionspace.getDim());
720    Data out(0.0,grad_shape,functionspace,true);    Data out(0.0,grad_shape,functionspace,true);
721    getDomain().setToGradient(out,*this);    getDomain()->setToGradient(out,*this);
722      blocktimer_increment("grad()", blocktimer_start);
723    return out;    return out;
724  }  }
725    
726  Data  Data
727  Data::grad() const  Data::grad() const
728  {  {
729    return gradOn(escript::function(getDomain()));    if (isEmpty())
730      {
731        throw DataException("Error - operation not permitted on instances of DataEmpty.");
732      }
733      return gradOn(escript::function(*getDomain()));
734  }  }
735    
736  int  int
737  Data::getDataPointSize() const  Data::getDataPointSize() const
738  {  {
739    return getPointDataView().noValues();    return m_data->getNoValues();
740  }  }
741    
742  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
743  Data::getLength() const  Data::getLength() const
744  {  {
745    return m_data->getLength();    return m_data->getLength();
746  }  }
747    
 const DataArrayView::ShapeType&  
 Data::getDataPointShape() const  
 {  
   return getPointDataView().getShape();  
 }  
   
748  const  const
749  boost::python::numeric::array  boost::python::numeric::array
750  Data::convertToNumArray()  Data:: getValueOfDataPoint(int dataPointNo)
751  {  {
752    //    size_t length=0;
753    // determine the total number of data points    int i, j, k, l;
   int numSamples = getNumSamples();  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
   int numDataPoints = numSamples * numDataPointsPerSample;  
   
754    //    //
755    // determine the rank and shape of each data point    // determine the rank and shape of each data point
756    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
757    DataArrayView::ShapeType dataPointShape = getDataPointShape();    const DataTypes::ShapeType& dataPointShape = getDataPointShape();
758    
759    //    //
760    // create the numeric array to be returned    // create the numeric array to be returned
761    boost::python::numeric::array numArray(0.0);    boost::python::numeric::array numArray(0.0);
762    
763    //    //
764    // the rank of the returned numeric array will be the rank of    // the shape of the returned numeric array will be the same
765    // the data points, plus one. Where the rank of the array is n,    // as that of the data point
766    // the last n-1 dimensions will be equal to the shape of the    int arrayRank = dataPointRank;
767    // data points, whilst the first dimension will be equal to the    const DataTypes::ShapeType& arrayShape = dataPointShape;
   // total number of data points. Thus the array will consist of  
   // a serial vector of the data points.  
   int arrayRank = dataPointRank + 1;  
   DataArrayView::ShapeType arrayShape;  
   arrayShape.push_back(numDataPoints);  
   for (int d=0; d<dataPointRank; d++) {  
      arrayShape.push_back(dataPointShape[d]);  
   }  
768    
769    //    //
770    // resize the numeric array to the shape just calculated    // resize the numeric array to the shape just calculated
771      if (arrayRank==0) {
772        numArray.resize(1);
773      }
774    if (arrayRank==1) {    if (arrayRank==1) {
775      numArray.resize(arrayShape[0]);      numArray.resize(arrayShape[0]);
776    }    }
# Line 490  Data::convertToNumArray() Line 783  Data::convertToNumArray()
783    if (arrayRank==4) {    if (arrayRank==4) {
784      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
785    }    }
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
   }  
786    
787    //    if (getNumDataPointsPerSample()>0) {
788    // loop through each data point in turn, loading the values for that data point         int sampleNo = dataPointNo/getNumDataPointsPerSample();
789    // into the numeric array.         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
790    int dataPoint = 0;         //
791    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {         // Check a valid sample number has been supplied
792      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
793        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
794        if (dataPointRank==0) {         }
795          numArray[dataPoint]=dataPointView();  
796        }         //
797        if (dataPointRank==1) {         // Check a valid data point number has been supplied
798          for (int i=0; i<dataPointShape[0]; i++) {         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
799            numArray[dataPoint][i]=dataPointView(i);             throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
800          }         }
801        }         // TODO: global error handling
802        if (dataPointRank==2) {         // create a view of the data if it is stored locally
803          for (int i=0; i<dataPointShape[0]; i++) {  //       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
804            for (int j=0; j<dataPointShape[1]; j++) {         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
805              numArray[dataPoint][i][j] = dataPointView(i,j);  
806            }  
807          }         switch( dataPointRank ){
808        }              case 0 :
809        if (dataPointRank==3) {                  numArray[0] = getDataAtOffset(offset);
810          for (int i=0; i<dataPointShape[0]; i++) {                  break;
811            for (int j=0; j<dataPointShape[1]; j++) {              case 1 :
812              for (int k=0; k<dataPointShape[2]; k++) {                  for( i=0; i<dataPointShape[0]; i++ )
813                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);                      numArray[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));
814              }                  break;
815            }              case 2 :
816          }                  for( i=0; i<dataPointShape[0]; i++ )
817        }                      for( j=0; j<dataPointShape[1]; j++)
818        if (dataPointRank==4) {                          numArray[make_tuple(i,j)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));
819          for (int i=0; i<dataPointShape[0]; i++) {                  break;
820            for (int j=0; j<dataPointShape[1]; j++) {              case 3 :
821              for (int k=0; k<dataPointShape[2]; k++) {                  for( i=0; i<dataPointShape[0]; i++ )
822                for (int l=0; l<dataPointShape[3]; l++) {                      for( j=0; j<dataPointShape[1]; j++ )
823                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);                          for( k=0; k<dataPointShape[2]; k++)
824                }                              numArray[make_tuple(i,j,k)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));
825              }                  break;
826            }              case 4 :
827          }                  for( i=0; i<dataPointShape[0]; i++ )
828        }                      for( j=0; j<dataPointShape[1]; j++ )
829        dataPoint++;                          for( k=0; k<dataPointShape[2]; k++ )
830      }                              for( l=0; l<dataPointShape[3]; l++)
831                                    numArray[make_tuple(i,j,k,l)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));
832                    break;
833        }
834    }    }
   
835    //    //
836    // return the loaded array    // return the array
837    return numArray;    return numArray;
838    
839  }  }
840    
841  const  void
842  boost::python::numeric::array  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
 Data::convertToNumArrayFromSampleNo(int sampleNo)  
843  {  {
844    //      // this will throw if the value cannot be represented
845    // Check a valid sample number has been supplied      boost::python::numeric::array num_array(py_object);
846    if (sampleNo >= getNumSamples()) {      setValueOfDataPointToArray(dataPointNo,num_array);
847      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  }
   }  
   
   //  
   // determine the number of data points per sample  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
   
   //  
   // determine the rank and shape of each data point  
   int dataPointRank = getDataPointRank();  
   DataArrayView::ShapeType dataPointShape = getDataPointShape();  
848    
849    void
850    Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)
851    {
852      if (isProtected()) {
853            throw DataException("Error - attempt to update protected Data object.");
854      }
855    //    //
856    // create the numeric array to be returned    // check rank
857    boost::python::numeric::array numArray(0.0);    if (num_array.getrank()<getDataPointRank())
858          throw DataException("Rank of numarray does not match Data object rank");
859    
860    //    //
861    // the rank of the returned numeric array will be the rank of    // check shape of num_array
862    // the data points, plus one. Where the rank of the array is n,    for (int i=0; i<getDataPointRank(); i++) {
863    // the last n-1 dimensions will be equal to the shape of the      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
864    // data points, whilst the first dimension will be equal to the         throw DataException("Shape of numarray does not match Data object rank");
   // total number of data points. Thus the array will consist of  
   // a serial vector of the data points.  
   int arrayRank = dataPointRank + 1;  
   DataArrayView::ShapeType arrayShape;  
   arrayShape.push_back(numDataPointsPerSample);  
   for (int d=0; d<dataPointRank; d++) {  
      arrayShape.push_back(dataPointShape[d]);  
865    }    }
   
866    //    //
867    // resize the numeric array to the shape just calculated    // make sure data is expanded:
868    if (arrayRank==1) {    //
869      numArray.resize(arrayShape[0]);    if (!isExpanded()) {
870    }      expand();
   if (arrayRank==2) {  
     numArray.resize(arrayShape[0],arrayShape[1]);  
   }  
   if (arrayRank==3) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);  
   }  
   if (arrayRank==4) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);  
871    }    }
872    if (arrayRank==5) {    if (getNumDataPointsPerSample()>0) {
873      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);         int sampleNo = dataPointNo/getNumDataPointsPerSample();
874           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
875           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
876      } else {
877           m_data->copyToDataPoint(-1, 0,num_array);
878    }    }
879    }
880    
881    //  void
882    // loop through each data point in turn, loading the values for that data point  Data::setValueOfDataPoint(int dataPointNo, const double value)
883    // into the numeric array.  {
884    for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {    if (isProtected()) {
885      DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);          throw DataException("Error - attempt to update protected Data object.");
     if (dataPointRank==0) {  
       numArray[dataPoint]=dataPointView();  
     }  
     if (dataPointRank==1) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         numArray[dataPoint][i]=dataPointView(i);  
       }  
     }  
     if (dataPointRank==2) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         for (int j=0; j<dataPointShape[1]; j++) {  
           numArray[dataPoint][i][j] = dataPointView(i,j);  
         }  
       }  
     }  
     if (dataPointRank==3) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         for (int j=0; j<dataPointShape[1]; j++) {  
           for (int k=0; k<dataPointShape[2]; k++) {  
             numArray[dataPoint][i][j][k]=dataPointView(i,j,k);  
           }  
         }  
       }  
     }  
     if (dataPointRank==4) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         for (int j=0; j<dataPointShape[1]; j++) {  
           for (int k=0; k<dataPointShape[2]; k++) {  
             for (int l=0; l<dataPointShape[3]; l++) {  
               numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);  
             }  
           }  
         }  
       }  
     }  
886    }    }
   
887    //    //
888    // return the loaded array    // make sure data is expanded:
889    return numArray;    if (!isExpanded()) {
890        expand();
891      }
892      if (getNumDataPointsPerSample()>0) {
893           int sampleNo = dataPointNo/getNumDataPointsPerSample();
894           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
895           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
896      } else {
897           m_data->copyToDataPoint(-1, 0,value);
898      }
899  }  }
900    
901  const  const
902  boost::python::numeric::array  boost::python::numeric::array
903  Data::convertToNumArrayFromDPNo(int sampleNo,  Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
                                 int dataPointNo)  
904  {  {
905    //    size_t length=0;
906    // Check a valid sample number has been supplied    int i, j, k, l, pos;
   if (sampleNo >= getNumSamples()) {  
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
   }  
   
   //  
   // Check a valid data point number has been supplied  
   if (dataPointNo >= getNumDataPointsPerSample()) {  
     throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");  
   }  
   
907    //    //
908    // determine the rank and shape of each data point    // determine the rank and shape of each data point
909    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
910    DataArrayView::ShapeType dataPointShape = getDataPointShape();    const DataTypes::ShapeType& dataPointShape = getDataPointShape();
911    
912    //    //
913    // create the numeric array to be returned    // create the numeric array to be returned
# Line 677  Data::convertToNumArrayFromDPNo(int samp Line 917  Data::convertToNumArrayFromDPNo(int samp
917    // the shape of the returned numeric array will be the same    // the shape of the returned numeric array will be the same
918    // as that of the data point    // as that of the data point
919    int arrayRank = dataPointRank;    int arrayRank = dataPointRank;
920    DataArrayView::ShapeType arrayShape = dataPointShape;    const DataTypes::ShapeType& arrayShape = dataPointShape;
921    
922    //    //
923    // resize the numeric array to the shape just calculated    // resize the numeric array to the shape just calculated
# Line 697  Data::convertToNumArrayFromDPNo(int samp Line 937  Data::convertToNumArrayFromDPNo(int samp
937      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
938    }    }
939    
940      // added for the MPI communication
941      length=1;
942      for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
943      double *tmpData = new double[length];
944    
945    //    //
946    // load the values for the data point into the numeric array.    // load the values for the data point into the numeric array.
   DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);  
   if (dataPointRank==0) {  
     numArray[0]=dataPointView();  
   }  
   if (dataPointRank==1) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       numArray[i]=dataPointView(i);  
     }  
   }  
   if (dataPointRank==2) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         numArray[i][j] = dataPointView(i,j);  
       }  
     }  
   }  
   if (dataPointRank==3) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         for (int k=0; k<dataPointShape[2]; k++) {  
           numArray[i][j][k]=dataPointView(i,j,k);  
         }  
       }  
     }  
   }  
   if (dataPointRank==4) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         for (int k=0; k<dataPointShape[2]; k++) {  
           for (int l=0; l<dataPointShape[3]; l++) {  
             numArray[i][j][k][l]=dataPointView(i,j,k,l);  
           }  
         }  
       }  
     }  
   }  
947    
948        // updated for the MPI case
949        if( get_MPIRank()==procNo ){
950                 if (getNumDataPointsPerSample()>0) {
951                    int sampleNo = dataPointNo/getNumDataPointsPerSample();
952                    int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
953                    //
954                    // Check a valid sample number has been supplied
955                    if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
956                      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
957                    }
958    
959                    //
960                    // Check a valid data point number has been supplied
961                    if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
962                      throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
963                    }
964                    // TODO: global error handling
965            // create a view of the data if it is stored locally
966            //DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
967            DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
968    
969            // pack the data from the view into tmpData for MPI communication
970            pos=0;
971            switch( dataPointRank ){
972                case 0 :
973                    tmpData[0] = getDataAtOffset(offset);
974                    break;
975                case 1 :
976                    for( i=0; i<dataPointShape[0]; i++ )
977                        tmpData[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));
978                    break;
979                case 2 :
980                    for( i=0; i<dataPointShape[0]; i++ )
981                        for( j=0; j<dataPointShape[1]; j++, pos++ )
982                            tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));
983                    break;
984                case 3 :
985                    for( i=0; i<dataPointShape[0]; i++ )
986                        for( j=0; j<dataPointShape[1]; j++ )
987                            for( k=0; k<dataPointShape[2]; k++, pos++ )
988                                tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));
989                    break;
990                case 4 :
991                    for( i=0; i<dataPointShape[0]; i++ )
992                        for( j=0; j<dataPointShape[1]; j++ )
993                            for( k=0; k<dataPointShape[2]; k++ )
994                                for( l=0; l<dataPointShape[3]; l++, pos++ )
995                                    tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));
996                    break;
997            }
998                }
999        }
1000            #ifdef PASO_MPI
1001            // broadcast the data to all other processes
1002        MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1003            #endif
1004    
1005        // unpack the data
1006        switch( dataPointRank ){
1007            case 0 :
1008                numArray[0]=tmpData[0];
1009                break;
1010            case 1 :
1011                for( i=0; i<dataPointShape[0]; i++ )
1012                    numArray[i]=tmpData[i];
1013                break;
1014            case 2 :
1015                for( i=0; i<dataPointShape[0]; i++ )
1016                    for( j=0; j<dataPointShape[1]; j++ )
1017                       numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
1018                break;
1019            case 3 :
1020                for( i=0; i<dataPointShape[0]; i++ )
1021                    for( j=0; j<dataPointShape[1]; j++ )
1022                        for( k=0; k<dataPointShape[2]; k++ )
1023                            numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
1024                break;
1025            case 4 :
1026                for( i=0; i<dataPointShape[0]; i++ )
1027                    for( j=0; j<dataPointShape[1]; j++ )
1028                        for( k=0; k<dataPointShape[2]; k++ )
1029                            for( l=0; l<dataPointShape[3]; l++ )
1030                                    numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
1031                break;
1032        }
1033    
1034        delete [] tmpData;
1035    //    //
1036    // return the loaded array    // return the loaded array
1037    return numArray;    return numArray;
1038  }  }
1039    
1040    
1041    
1042  boost::python::numeric::array  boost::python::numeric::array
1043  Data::integrate() const  Data::integrate() const
1044  {  {
1045    int index;    int index;
1046    int rank = getDataPointRank();    int rank = getDataPointRank();
1047    DataArrayView::ShapeType shape = getDataPointShape();    DataTypes::ShapeType shape = getDataPointShape();
1048      int dataPointSize = getDataPointSize();
1049    
1050    //    //
1051    // calculate the integral values    // calculate the integral values
1052    vector<double> integrals(getDataPointSize());    vector<double> integrals(dataPointSize);
1053    AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);    vector<double> integrals_local(dataPointSize);
1054    #ifdef PASO_MPI
1055      AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals_local,*this);
1056      // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1057      double *tmp = new double[dataPointSize];
1058      double *tmp_local = new double[dataPointSize];
1059      for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1060      MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1061      for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1062      delete[] tmp;
1063      delete[] tmp_local;
1064    #else
1065      AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);
1066    #endif
1067    
1068    //    //
1069    // create the numeric array to be returned    // create the numeric array to be returned
# Line 770  Data::integrate() const Line 1082  Data::integrate() const
1082      }      }
1083    }    }
1084    if (rank==2) {    if (rank==2) {
1085      bp_array.resize(shape[0],shape[1]);         bp_array.resize(shape[0],shape[1]);
1086      for (int i=0; i<shape[0]; i++) {         for (int i=0; i<shape[0]; i++) {
1087        for (int j=0; j<shape[1]; j++) {           for (int j=0; j<shape[1]; j++) {
1088          index = i + shape[0] * j;             index = i + shape[0] * j;
1089          bp_array[i,j] = integrals[index];             bp_array[make_tuple(i,j)] = integrals[index];
1090        }           }
1091      }         }
1092    }    }
1093    if (rank==3) {    if (rank==3) {
1094      bp_array.resize(shape[0],shape[1],shape[2]);      bp_array.resize(shape[0],shape[1],shape[2]);
# Line 784  Data::integrate() const Line 1096  Data::integrate() const
1096        for (int j=0; j<shape[1]; j++) {        for (int j=0; j<shape[1]; j++) {
1097          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
1098            index = i + shape[0] * ( j + shape[1] * k );            index = i + shape[0] * ( j + shape[1] * k );
1099            bp_array[i,j,k] = integrals[index];            bp_array[make_tuple(i,j,k)] = integrals[index];
1100          }          }
1101        }        }
1102      }      }
# Line 796  Data::integrate() const Line 1108  Data::integrate() const
1108          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
1109            for (int l=0; l<shape[3]; l++) {            for (int l=0; l<shape[3]; l++) {
1110              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
1111              bp_array[i,j,k,l] = integrals[index];              bp_array[make_tuple(i,j,k,l)] = integrals[index];
1112            }            }
1113          }          }
1114        }        }
# Line 811  Data::integrate() const Line 1123  Data::integrate() const
1123  Data  Data
1124  Data::sin() const  Data::sin() const
1125  {  {
1126    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1127  }  }
1128    
1129  Data  Data
1130  Data::cos() const  Data::cos() const
1131  {  {
1132    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1133  }  }
1134    
1135  Data  Data
1136  Data::tan() const  Data::tan() const
1137  {  {
1138    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1139  }  }
1140    
1141  Data  Data
1142  Data::log() const  Data::asin() const
1143    {
1144      return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1145    }
1146    
1147    Data
1148    Data::acos() const
1149    {
1150      return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1151    }
1152    
1153    
1154    Data
1155    Data::atan() const
1156    {
1157      return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1158    }
1159    
1160    Data
1161    Data::sinh() const
1162    {
1163        return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1164    
1165    }
1166    
1167    Data
1168    Data::cosh() const
1169    {
1170        return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1171    }
1172    
1173    Data
1174    Data::tanh() const
1175    {
1176        return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1177    }
1178    
1179    
1180    Data
1181    Data::erf() const
1182    {
1183    #ifdef _WIN32
1184      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1185    #else
1186      return C_TensorUnaryOperation(*this, ::erf);
1187    #endif
1188    }
1189    
1190    Data
1191    Data::asinh() const
1192    {
1193    #ifdef _WIN32
1194      return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1195    #else
1196      return C_TensorUnaryOperation(*this, ::asinh);
1197    #endif
1198    }
1199    
1200    Data
1201    Data::acosh() const
1202    {
1203    #ifdef _WIN32
1204      return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1205    #else
1206      return C_TensorUnaryOperation(*this, ::acosh);
1207    #endif
1208    }
1209    
1210    Data
1211    Data::atanh() const
1212  {  {
1213    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);  #ifdef _WIN32
1214      return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1215    #else
1216      return C_TensorUnaryOperation(*this, ::atanh);
1217    #endif
1218  }  }
1219    
1220  Data  Data
1221  Data::ln() const  Data::log10() const
1222  {  {
1223    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1224    }
1225    
1226    Data
1227    Data::log() const
1228    {
1229      return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1230  }  }
1231    
1232  Data  Data
1233  Data::sign() const  Data::sign() const
1234  {  {
1235    return escript::unaryOp(*this,escript::fsign);    return C_TensorUnaryOperation(*this, escript::fsign);
1236  }  }
1237    
1238  Data  Data
1239  Data::abs() const  Data::abs() const
1240  {  {
1241    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1242  }  }
1243    
1244  Data  Data
1245  Data::neg() const  Data::neg() const
1246  {  {
1247    return escript::unaryOp(*this,negate<double>());    return C_TensorUnaryOperation(*this, negate<double>());
1248  }  }
1249    
1250  Data  Data
1251  Data::pos() const  Data::pos() const
1252  {  {
1253    return (*this);    Data result;
1254      // perform a deep copy
1255      result.copy(*this);
1256      return result;
1257  }  }
1258    
1259  Data  Data
1260  Data::exp() const  Data::exp() const
1261  {  {
1262    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1263  }  }
1264    
1265  Data  Data
1266  Data::sqrt() const  Data::sqrt() const
1267  {  {
1268    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1269  }  }
1270    
1271  double  double
1272  Data::Lsup() const  Data::Lsup() const
1273  {  {
1274      double localValue;
1275    //    //
1276    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
   return algorithm(DataAlgorithmAdapter<AbsMax>(0));  
 }  
1277    
1278  double    AbsMax abs_max_func;
1279  Data::Linf() const    localValue = algorithm(abs_max_func,0);
1280  {  #ifdef PASO_MPI
1281    //    double globalValue;
1282    // set the initial absolute minimum value to max double    MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1283    return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));    return globalValue;
1284    #else
1285      return localValue;
1286    #endif
1287  }  }
1288    
1289  double  double
1290  Data::sup() const  Data::sup() const
1291  {  {
1292      double localValue;
1293    //    //
1294    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1295    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1296      localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1297    #ifdef PASO_MPI
1298      double globalValue;
1299      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1300      return globalValue;
1301    #else
1302      return localValue;
1303    #endif
1304  }  }
1305    
1306  double  double
1307  Data::inf() const  Data::inf() const
1308  {  {
1309      double localValue;
1310    //    //
1311    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1312    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1313      localValue = algorithm(fmin_func,numeric_limits<double>::max());
1314    #ifdef PASO_MPI
1315      double globalValue;
1316      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1317      return globalValue;
1318    #else
1319      return localValue;
1320    #endif
1321  }  }
1322    
1323    /* TODO */
1324    /* global reduction */
1325  Data  Data
1326  Data::maxval() const  Data::maxval() const
1327  {  {
1328    //    //
1329    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1330    return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1331      return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1332  }  }
1333    
1334  Data  Data
# Line 919  Data::minval() const Line 1336  Data::minval() const
1336  {  {
1337    //    //
1338    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1339    return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1340      return dp_algorithm(fmin_func,numeric_limits<double>::max());
1341  }  }
1342    
1343  const boost::python::tuple  Data
1344  Data::mindp() const  Data::swapaxes(const int axis0, const int axis1) const
1345  {  {
1346    Data temp=minval();       int axis0_tmp,axis1_tmp;
1347         DataTypes::ShapeType s=getDataPointShape();
1348    int numSamples=temp.getNumSamples();       DataTypes::ShapeType ev_shape;
1349    int numDPPSample=temp.getNumDataPointsPerSample();       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1350         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1351    int i,j,lowi=0,lowj=0;       int rank=getDataPointRank();
1352    double min=numeric_limits<double>::max();       if (rank<2) {
1353            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1354    for (i=0; i<numSamples; i++) {       }
1355      for (j=0; j<numDPPSample; j++) {       if (axis0<0 || axis0>rank-1) {
1356        double next=temp.getDataPoint(i,j)();          throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1357        if (next<min) {       }
1358          min=next;       if (axis1<0 || axis1>rank-1) {
1359          lowi=i;           throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1360          lowj=j;       }
1361        }       if (axis0 == axis1) {
1362      }           throw DataException("Error - Data::swapaxes: axis indices must be different.");
1363    }       }
1364         if (axis0 > axis1) {
1365             axis0_tmp=axis1;
1366             axis1_tmp=axis0;
1367         } else {
1368             axis0_tmp=axis0;
1369             axis1_tmp=axis1;
1370         }
1371         for (int i=0; i<rank; i++) {
1372           if (i == axis0_tmp) {
1373              ev_shape.push_back(s[axis1_tmp]);
1374           } else if (i == axis1_tmp) {
1375              ev_shape.push_back(s[axis0_tmp]);
1376           } else {
1377              ev_shape.push_back(s[i]);
1378           }
1379         }
1380         Data ev(0.,ev_shape,getFunctionSpace());
1381         ev.typeMatchRight(*this);
1382         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1383         return ev;
1384    
1385    }
1386    
1387    Data
1388    Data::symmetric() const
1389    {
1390         // check input
1391         DataTypes::ShapeType s=getDataPointShape();
1392         if (getDataPointRank()==2) {
1393            if(s[0] != s[1])
1394               throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1395         }
1396         else if (getDataPointRank()==4) {
1397            if(!(s[0] == s[2] && s[1] == s[3]))
1398               throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1399         }
1400         else {
1401            throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1402         }
1403         Data ev(0.,getDataPointShape(),getFunctionSpace());
1404         ev.typeMatchRight(*this);
1405         m_data->symmetric(ev.m_data.get());
1406         return ev;
1407    }
1408    
1409    Data
1410    Data::nonsymmetric() const
1411    {
1412         // check input
1413         DataTypes::ShapeType s=getDataPointShape();
1414         if (getDataPointRank()==2) {
1415            if(s[0] != s[1])
1416               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1417            DataTypes::ShapeType ev_shape;
1418            ev_shape.push_back(s[0]);
1419            ev_shape.push_back(s[1]);
1420            Data ev(0.,ev_shape,getFunctionSpace());
1421            ev.typeMatchRight(*this);
1422            m_data->nonsymmetric(ev.m_data.get());
1423            return ev;
1424         }
1425         else if (getDataPointRank()==4) {
1426            if(!(s[0] == s[2] && s[1] == s[3]))
1427               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1428            DataTypes::ShapeType ev_shape;
1429            ev_shape.push_back(s[0]);
1430            ev_shape.push_back(s[1]);
1431            ev_shape.push_back(s[2]);
1432            ev_shape.push_back(s[3]);
1433            Data ev(0.,ev_shape,getFunctionSpace());
1434            ev.typeMatchRight(*this);
1435            m_data->nonsymmetric(ev.m_data.get());
1436            return ev;
1437         }
1438         else {
1439            throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1440         }
1441    }
1442    
1443    Data
1444    Data::trace(int axis_offset) const
1445    {
1446         DataTypes::ShapeType s=getDataPointShape();
1447         if (getDataPointRank()==2) {
1448            DataTypes::ShapeType ev_shape;
1449            Data ev(0.,ev_shape,getFunctionSpace());
1450            ev.typeMatchRight(*this);
1451            m_data->trace(ev.m_data.get(), axis_offset);
1452            return ev;
1453         }
1454         if (getDataPointRank()==3) {
1455            DataTypes::ShapeType ev_shape;
1456            if (axis_offset==0) {
1457              int s2=s[2];
1458              ev_shape.push_back(s2);
1459            }
1460            else if (axis_offset==1) {
1461              int s0=s[0];
1462              ev_shape.push_back(s0);
1463            }
1464            Data ev(0.,ev_shape,getFunctionSpace());
1465            ev.typeMatchRight(*this);
1466            m_data->trace(ev.m_data.get(), axis_offset);
1467            return ev;
1468         }
1469         if (getDataPointRank()==4) {
1470            DataTypes::ShapeType ev_shape;
1471            if (axis_offset==0) {
1472              ev_shape.push_back(s[2]);
1473              ev_shape.push_back(s[3]);
1474            }
1475            else if (axis_offset==1) {
1476              ev_shape.push_back(s[0]);
1477              ev_shape.push_back(s[3]);
1478            }
1479        else if (axis_offset==2) {
1480          ev_shape.push_back(s[0]);
1481          ev_shape.push_back(s[1]);
1482        }
1483            Data ev(0.,ev_shape,getFunctionSpace());
1484            ev.typeMatchRight(*this);
1485        m_data->trace(ev.m_data.get(), axis_offset);
1486            return ev;
1487         }
1488         else {
1489            throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1490         }
1491    }
1492    
1493    return make_tuple(lowi,lowj);  Data
1494    Data::transpose(int axis_offset) const
1495    {
1496         DataTypes::ShapeType s=getDataPointShape();
1497         DataTypes::ShapeType ev_shape;
1498         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1499         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1500         int rank=getDataPointRank();
1501         if (axis_offset<0 || axis_offset>rank) {
1502            throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1503         }
1504         for (int i=0; i<rank; i++) {
1505           int index = (axis_offset+i)%rank;
1506           ev_shape.push_back(s[index]); // Append to new shape
1507         }
1508         Data ev(0.,ev_shape,getFunctionSpace());
1509         ev.typeMatchRight(*this);
1510         m_data->transpose(ev.m_data.get(), axis_offset);
1511         return ev;
1512  }  }
1513    
1514  Data  Data
1515  Data::length() const  Data::eigenvalues() const
1516    {
1517         // check input
1518         DataTypes::ShapeType s=getDataPointShape();
1519         if (getDataPointRank()!=2)
1520            throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1521         if(s[0] != s[1])
1522            throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1523         // create return
1524         DataTypes::ShapeType ev_shape(1,s[0]);
1525         Data ev(0.,ev_shape,getFunctionSpace());
1526         ev.typeMatchRight(*this);
1527         m_data->eigenvalues(ev.m_data.get());
1528         return ev;
1529    }
1530    
1531    const boost::python::tuple
1532    Data::eigenvalues_and_eigenvectors(const double tol) const
1533  {  {
1534    return dp_algorithm(DataAlgorithmAdapter<Length>(0));       DataTypes::ShapeType s=getDataPointShape();
1535         if (getDataPointRank()!=2)
1536            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1537         if(s[0] != s[1])
1538            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1539         // create return
1540         DataTypes::ShapeType ev_shape(1,s[0]);
1541         Data ev(0.,ev_shape,getFunctionSpace());
1542         ev.typeMatchRight(*this);
1543         DataTypes::ShapeType V_shape(2,s[0]);
1544         Data V(0.,V_shape,getFunctionSpace());
1545         V.typeMatchRight(*this);
1546         m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1547         return make_tuple(boost::python::object(ev),boost::python::object(V));
1548  }  }
1549    
1550  Data  const boost::python::tuple
1551  Data::trace() const  Data::minGlobalDataPoint() const
1552  {  {
1553    return dp_algorithm(DataAlgorithmAdapter<Trace>(0));    // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1554      // abort (for unknown reasons) if there are openmp directives with it in the
1555      // surrounding function
1556    
1557      int DataPointNo;
1558      int ProcNo;
1559      calc_minGlobalDataPoint(ProcNo,DataPointNo);
1560      return make_tuple(ProcNo,DataPointNo);
1561  }  }
1562    
1563  Data  void
1564  Data::transpose(int axis) const  Data::calc_minGlobalDataPoint(int& ProcNo,
1565                            int& DataPointNo) const
1566  {  {
1567    // not implemented    int i,j;
1568    throw DataException("Error - Data::transpose not implemented yet.");    int lowi=0,lowj=0;
1569    return Data();    double min=numeric_limits<double>::max();
1570    
1571      Data temp=minval();
1572    
1573      int numSamples=temp.getNumSamples();
1574      int numDPPSample=temp.getNumDataPointsPerSample();
1575    
1576      double next,local_min;
1577      int local_lowi,local_lowj;
1578    
1579      #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1580      {
1581        local_min=min;
1582        #pragma omp for private(i,j) schedule(static)
1583        for (i=0; i<numSamples; i++) {
1584          for (j=0; j<numDPPSample; j++) {
1585            next=temp.getDataAtOffset(temp.getDataOffset(i,j));
1586            if (next<local_min) {
1587              local_min=next;
1588              local_lowi=i;
1589              local_lowj=j;
1590            }
1591          }
1592        }
1593        #pragma omp critical
1594        if (local_min<min) {
1595          min=local_min;
1596          lowi=local_lowi;
1597          lowj=local_lowj;
1598        }
1599      }
1600    
1601    #ifdef PASO_MPI
1602        // determine the processor on which the minimum occurs
1603        next = temp.getDataPoint(lowi,lowj);
1604        int lowProc = 0;
1605        double *globalMins = new double[get_MPISize()+1];
1606        int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1607    
1608        if( get_MPIRank()==0 ){
1609            next = globalMins[lowProc];
1610            for( i=1; i<get_MPISize(); i++ )
1611                if( next>globalMins[i] ){
1612                    lowProc = i;
1613                    next = globalMins[i];
1614                }
1615        }
1616        MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1617    
1618        delete [] globalMins;
1619        ProcNo = lowProc;
1620    #else
1621        ProcNo = 0;
1622    #endif
1623      DataPointNo = lowj + lowi * numDPPSample;
1624  }  }
1625    
1626  void  void
1627  Data::saveDX(std::string fileName) const  Data::saveDX(std::string fileName) const
1628  {  {
1629    getDomain().saveDX(fileName,*this);    if (isEmpty())
1630      {
1631        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1632      }
1633      boost::python::dict args;
1634      args["data"]=boost::python::object(this);
1635      getDomain()->saveDX(fileName,args);
1636    return;    return;
1637  }  }
1638    
1639  void  void
1640  Data::saveVTK(std::string fileName) const  Data::saveVTK(std::string fileName) const
1641  {  {
1642    getDomain().saveVTK(fileName,*this);    if (isEmpty())
1643      {
1644        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1645      }
1646      boost::python::dict args;
1647      args["data"]=boost::python::object(this);
1648      getDomain()->saveVTK(fileName,args);
1649    return;    return;
1650  }  }
1651    
1652  Data&  Data&
1653  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1654  {  {
1655      if (isProtected()) {
1656            throw DataException("Error - attempt to update protected Data object.");
1657      }
1658    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1659    return (*this);    return (*this);
1660  }  }
# Line 991  Data::operator+=(const Data& right) Line 1662  Data::operator+=(const Data& right)
1662  Data&  Data&
1663  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1664  {  {
1665    binaryOp(right,plus<double>());    Data tmp(right,getFunctionSpace(),false);
1666      binaryOp(tmp,plus<double>());
1667      return (*this);
1668    }
1669    Data&
1670    Data::operator=(const Data& other)
1671    {
1672      copy(other);
1673    return (*this);    return (*this);
1674  }  }
1675    
1676  Data&  Data&
1677  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1678  {  {
1679      if (isProtected()) {
1680            throw DataException("Error - attempt to update protected Data object.");
1681      }
1682    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1683    return (*this);    return (*this);
1684  }  }
# Line 1005  Data::operator-=(const Data& right) Line 1686  Data::operator-=(const Data& right)
1686  Data&  Data&
1687  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1688  {  {
1689    binaryOp(right,minus<double>());    Data tmp(right,getFunctionSpace(),false);
1690      binaryOp(tmp,minus<double>());
1691    return (*this);    return (*this);
1692  }  }
1693    
1694  Data&  Data&
1695  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1696  {  {
1697      if (isProtected()) {
1698            throw DataException("Error - attempt to update protected Data object.");
1699      }
1700    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1701    return (*this);    return (*this);
1702  }  }
# Line 1019  Data::operator*=(const Data& right) Line 1704  Data::operator*=(const Data& right)
1704  Data&  Data&
1705  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1706  {  {
1707    binaryOp(right,multiplies<double>());    Data tmp(right,getFunctionSpace(),false);
1708      binaryOp(tmp,multiplies<double>());
1709    return (*this);    return (*this);
1710  }  }
1711    
1712  Data&  Data&
1713  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1714  {  {
1715      if (isProtected()) {
1716            throw DataException("Error - attempt to update protected Data object.");
1717      }
1718    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1719    return (*this);    return (*this);
1720  }  }
# Line 1033  Data::operator/=(const Data& right) Line 1722  Data::operator/=(const Data& right)
1722  Data&  Data&
1723  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1724  {  {
1725    binaryOp(right,divides<double>());    Data tmp(right,getFunctionSpace(),false);
1726      binaryOp(tmp,divides<double>());
1727    return (*this);    return (*this);
1728  }  }
1729    
1730  Data  Data
1731    Data::rpowO(const boost::python::object& left) const
1732    {
1733      Data left_d(left,*this);
1734      return left_d.powD(*this);
1735    }
1736    
1737    Data
1738  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1739  {  {
1740    Data result;    Data tmp(right,getFunctionSpace(),false);
1741    result.copy(*this);    return powD(tmp);
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
1742  }  }
1743    
1744  Data  Data
1745  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1746  {  {
1747    Data result;    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
   result.copy(*this);  
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
1748  }  }
1749    
1750  //  //
1751  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1752  Data  Data
1753  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
1754  {  {
1755    Data result;    return C_TensorBinaryOperation(left, right, plus<double>());
   //  
   // perform a deep copy  
   result.copy(left);  
   result+=right;  
   return result;  
1756  }  }
1757    
1758  //  //
1759  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1760  Data  Data
1761  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
1762  {  {
1763    Data result;    return C_TensorBinaryOperation(left, right, minus<double>());
   //  
   // perform a deep copy  
   result.copy(left);  
   result-=right;  
   return result;  
1764  }  }
1765    
1766  //  //
1767  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1768  Data  Data
1769  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
1770  {  {
1771    Data result;    return C_TensorBinaryOperation(left, right, multiplies<double>());
   //  
   // perform a deep copy  
   result.copy(left);  
   result*=right;  
   return result;  
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  {  {
1779    Data result;    return C_TensorBinaryOperation(left, right, divides<double>());
   //  
   // perform a deep copy  
   result.copy(left);  
   result/=right;  
   return result;  
1780  }  }
1781    
1782  //  //
1783  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1784  Data  Data
1785  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1786  {  {
1787    //    return left+Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result+=right;  
   return result;  
1788  }  }
1789    
1790  //  //
1791  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1792  Data  Data
1793  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1794  {  {
1795    //    return left-Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result-=right;  
   return result;  
1796  }  }
1797    
1798  //  //
1799  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1800  Data  Data
1801  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1802  {  {
1803    //    return left*Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result*=right;  
   return result;  
1804  }  }
1805    
1806  //  //
1807  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1808  Data  Data
1809  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1810  {  {
1811    //    return left/Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result/=right;  
   return result;  
1812  }  }
1813    
1814  //  //
1815  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1816  Data  Data
1817  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1818  {  {
1819    //    return Data(left,right.getFunctionSpace(),false)+right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result+=right;  
   return result;  
1820  }  }
1821    
1822  //  //
1823  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1824  Data  Data
1825  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1826  {  {
1827    //    return Data(left,right.getFunctionSpace(),false)-right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result-=right;  
   return result;  
1828  }  }
1829    
1830  //  //
1831  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1832  Data  Data
1833  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1834  {  {
1835    //    return Data(left,right.getFunctionSpace(),false)*right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result*=right;  
   return result;  
1836  }  }
1837    
1838  //  //
1839  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1840  Data  Data
1841  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1842  {  {
1843    //    return Data(left,right.getFunctionSpace(),false)/right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result/=right;  
   return result;  
1844  }  }
1845    
1846  //  //
 // NOTE: It is essential to specify the namepsace this operator belongs to  
1847  //bool escript::operator==(const Data& left, const Data& right)  //bool escript::operator==(const Data& left, const Data& right)
1848  //{  //{
1849  //  /*  //  /*
# Line 1269  escript::operator/(const boost::python:: Line 1888  escript::operator/(const boost::python::
1888  //  return ret;  //  return ret;
1889  //}  //}
1890    
1891    /* TODO */
1892    /* global reduction */
1893  Data  Data
1894  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
1895  {  {
1896    const DataArrayView& view=getPointDataView();  //  const DataArrayView& view=getPointDataView();
1897    
1898    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
1899    
1900    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=getDataPointRank()) {
1901      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
1902    }    }
1903    
1904    return getSlice(slice_region);    return getSlice(slice_region);
1905  }  }
1906    
1907    /* TODO */
1908    /* global reduction */
1909  Data  Data
1910  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataTypes::RegionType& region) const
1911  {  {
1912    return Data(*this,region);    return Data(*this,region);
1913  }  }
1914    
1915    /* TODO */
1916    /* global reduction */
1917  void  void
1918  Data::setItemO(const boost::python::object& key,  Data::setItemO(const boost::python::object& key,
1919                 const boost::python::object& value)                 const boost::python::object& value)
# Line 1301  void Line 1926  void
1926  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
1927                 const Data& value)                 const Data& value)
1928  {  {
1929    const DataArrayView& view=getPointDataView();  //  const DataArrayView& view=getPointDataView();
1930    DataArrayView::RegionType slice_region=view.getSliceRegion(key);  
1931    if (slice_region.size()!=view.getRank()) {    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
1932      if (slice_region.size()!=getDataPointRank()) {
1933      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
1934    }    }
1935    if (getFunctionSpace()!=value.getFunctionSpace()) {    if (getFunctionSpace()!=value.getFunctionSpace()) {
# Line 1315  Data::setItemD(const boost::python::obje Line 1941  Data::setItemD(const boost::python::obje
1941    
1942  void  void
1943  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1944                 const DataArrayView::RegionType& region)                 const DataTypes::RegionType& region)
1945  {  {
1946      if (isProtected()) {
1947            throw DataException("Error - attempt to update protected Data object.");
1948      }
1949    Data tempValue(value);    Data tempValue(value);
1950    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
1951    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 1352  Data::typeMatchRight(const Data& right) Line 1981  Data::typeMatchRight(const Data& right)
1981  }  }
1982    
1983  void  void
1984    Data::setTaggedValueByName(std::string name,
1985                               const boost::python::object& value)
1986    {
1987         if (getFunctionSpace().getDomain()->isValidTagName(name)) {
1988            int tagKey=getFunctionSpace().getDomain()->getTag(name);
1989            setTaggedValue(tagKey,value);
1990         }
1991    }
1992    void
1993  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
1994                       const boost::python::object& value)                       const boost::python::object& value)
1995  {  {
1996      if (isProtected()) {
1997            throw DataException("Error - attempt to update protected Data object.");
1998      }
1999    //    //
2000    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2001    tag();    if (isConstant()) tag();
2002    
2003    if (!isTagged()) {    numeric::array asNumArray(value);
2004      throw DataException("Error - DataTagged conversion failed!!");  
2005    
2006      // extract the shape of the numarray
2007      DataTypes::ShapeType tempShape;
2008      for (int i=0; i < asNumArray.getrank(); i++) {
2009        tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
2010    }    }
2011    
2012    //    // get the space for the data vector
2013    // Construct DataArray from boost::python::object input value  //   int len = DataTypes::noValues(tempShape);
2014    DataArray valueDataArray(value);  //   DataVector temp_data(len, 0.0, len);
2015    //   DataArrayView temp_dataView(temp_data, tempShape);
2016    //   temp_dataView.copy(asNumArray);
2017    
2018      DataVector temp_data2;
2019      temp_data2.copyFromNumArray(asNumArray);
2020    
2021    //    //
2022    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
2023    m_data->setTaggedValue(tagKey,valueDataArray.getView());    //m_data->setTaggedValue(tagKey,temp_dataView);
2024    
2025        m_data->setTaggedValue(tagKey,tempShape, temp_data2);
2026  }  }
2027    
2028    
2029  void  void
2030  Data::setTaggedValueFromCPP(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
2031                              const DataArrayView& value)                  const DataTypes::ShapeType& pointshape,
2032                                const DataTypes::ValueType& value,
2033                    int dataOffset)
2034  {  {
2035      if (isProtected()) {
2036            throw DataException("Error - attempt to update protected Data object.");
2037      }
2038    //    //
2039    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2040    tag();    if (isConstant()) tag();
2041    
   if (!isTagged()) {  
     throw DataException("Error - DataTagged conversion failed!!");  
   }  
                                                                                                                 
2042    //    //
2043    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
2044    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
2045  }  }
2046    
2047  void  int
2048  Data::setRefValue(int ref,  Data::getTagNumber(int dpno)
                   const boost::python::numeric::array& value)  
2049  {  {
2050    //    if (isEmpty())
2051    // Construct DataArray from boost::python::object input value    {
2052    DataArray valueDataArray(value);      throw DataException("Error - operation not permitted on instances of DataEmpty.");
2053      }
2054    //    return getFunctionSpace().getTagFromDataPointNo(dpno);
   // Call DataAbstract::setRefValue  
   m_data->setRefValue(ref,valueDataArray);  
2055  }  }
2056    
 void  
 Data::getRefValue(int ref,  
                   boost::python::numeric::array& value)  
 {  
   //  
   // Construct DataArray for boost::python::object return value  
   DataArray valueDataArray(value);  
   
   //  
   // Load DataArray with values from data-points specified by ref  
   m_data->getRefValue(ref,valueDataArray);  
2057    
2058    //  ostream& escript::operator<<(ostream& o, const Data& data)
2059    // Load values from valueDataArray into return numarray  {
2060      o << data.toString();
2061    // extract the shape of the numarray    return o;
2062    int rank = value.getrank();  }
   DataArrayView::ShapeType shape;  
   for (int i=0; i < rank; i++) {  
     shape.push_back(extract<int>(value.getshape()[i]));  
   }  
   
   // and load the numarray with the data from the DataArray  
   DataArrayView valueView = valueDataArray.getView();  
2063    
2064    if (rank==0) {  Data
2065      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");  escript::C_GeneralTensorProduct(Data& arg_0,
2066    }                       Data& arg_1,
2067    if (rank==1) {                       int axis_offset,
2068      for (int i=0; i < shape[0]; i++) {                       int transpose)
2069        value[i] = valueView(i);  {
2070      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2071      // SM is the product of the last axis_offset entries in arg_0.getShape().
2072    
2073      // Interpolate if necessary and find an appropriate function space
2074      Data arg_0_Z, arg_1_Z;
2075      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2076        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2077          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2078          arg_1_Z = Data(arg_1);
2079        }
2080        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2081          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2082          arg_0_Z =Data(arg_0);
2083      }      }
2084        else {
2085          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2086        }
2087      } else {
2088          arg_0_Z = Data(arg_0);
2089          arg_1_Z = Data(arg_1);
2090    }    }
2091    if (rank==2) {    // Get rank and shape of inputs
2092      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");    int rank0 = arg_0_Z.getDataPointRank();
2093    }    int rank1 = arg_1_Z.getDataPointRank();
2094    if (rank==3) {    const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2095      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");    const DataTypes::ShapeType& shape1 = arg_1_Z.getDataPointShape();
2096    }  
2097    if (rank==4) {    // Prepare for the loops of the product and verify compatibility of shapes
2098      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");    int start0=0, start1=0;
2099    }    if (transpose == 0)       {}
2100      else if (transpose == 1)  { start0 = axis_offset; }
2101  }    else if (transpose == 2)  { start1 = rank1-axis_offset; }
2102      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2103  void  
2104  Data::archiveData(const std::string fileName)  
2105  {    // Adjust the shapes for transpose
2106    cout << "Archiving Data object to: " << fileName << endl;    DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather
2107      DataTypes::ShapeType tmpShape1(rank1);    // than using push_back
2108    //    for (int i=0; i<rank0; i++)   { tmpShape0[i]=shape0[(i+start0)%rank0]; }
2109    // Determine type of this Data object    for (int i=0; i<rank1; i++)   { tmpShape1[i]=shape1[(i+start1)%rank1]; }
2110    int dataType = -1;  
2111    #if 0
2112      // For debugging: show shape after transpose
2113      char tmp[100];
2114      std::string shapeStr;
2115      shapeStr = "(";
2116      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2117      shapeStr += ")";
2118      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2119      shapeStr = "(";
2120      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2121      shapeStr += ")";
2122      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2123    #endif
2124    
2125      // Prepare for the loops of the product
2126      int SL=1, SM=1, SR=1;
2127      for (int i=0; i<rank0-axis_offset; i++)   {
2128        SL *= tmpShape0[i];
2129      }
2130      for (int i=rank0-axis_offset; i<rank0; i++)   {
2131        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2132          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2133        }
2134        SM *= tmpShape0[i];
2135      }
2136      for (int i=axis_offset; i<rank1; i++)     {
2137        SR *= tmpShape1[i];
2138      }
2139    
2140      // Define the shape of the output (rank of shape is the sum of the loop ranges below)
2141      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);  
2142      {         // block to limit the scope of out_index
2143         int out_index=0;
2144         for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
2145         for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2146      }
2147    
2148      // Declare output Data object
2149      Data res;
2150    
2151      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2152        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2153        double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));
2154        double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));
2155        double *ptr_2 = &(res.getDataAtOffset(0));
2156        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2157      }
2158      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2159    
2160        // Prepare the DataConstant input
2161        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2162        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2163    
2164        // Borrow DataTagged input from Data object
2165        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2166        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2167    
2168        // Prepare a DataTagged output 2
2169        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2170        res.tag();
2171        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2172        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2173    
2174        // Prepare offset into DataConstant
2175        int offset_0 = tmp_0->getPointOffset(0,0);
2176        double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2177        // Get the views
2178    //     DataArrayView view_1 = tmp_1->getDefaultValue();
2179    //     DataArrayView view_2 = tmp_2->getDefaultValue();
2180    //     // Get the pointers to the actual data
2181    //     double *ptr_1 = &((view_1.getData())[0]);
2182    //     double *ptr_2 = &((view_2.getData())[0]);
2183    
2184        double *ptr_1 = &(tmp_1->getDefaultValue(0));
2185        double *ptr_2 = &(tmp_2->getDefaultValue(0));
2186    
2187    
2188        // Compute an MVP for the default
2189        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2190        // Compute an MVP for each tag
2191        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2192        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2193        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2194          tmp_2->addTag(i->first);
2195    //       DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2196    //       DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2197    //       double *ptr_1 = &view_1.getData(0);
2198    //       double *ptr_2 = &view_2.getData(0);
2199    
2200          double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2201          double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2202        
2203          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2204        }
2205    
2206      }
2207      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2208    
2209        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2210        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2211        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2212        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2213        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2214        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2215        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2216        int sampleNo_1,dataPointNo_1;
2217        int numSamples_1 = arg_1_Z.getNumSamples();
2218        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2219        int offset_0 = tmp_0->getPointOffset(0,0);
2220        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2221        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2222          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2223            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2224            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2225            double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2226            double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2227            double *ptr_2 = &(res.getDataAtOffset(offset_2));
2228            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2229          }
2230        }
2231    
   if (isEmpty()) {  
     dataType = 0;  
     cout << "\tdataType: DataEmpty" << endl;  
   }  
   if (isConstant()) {  
     dataType = 1;  
     cout << "\tdataType: DataConstant" << endl;  
   }  
   if (isTagged()) {  
     dataType = 2;  
     cout << "\tdataType: DataTagged" << endl;  
   }  
   if (isExpanded()) {  
     dataType = 3;  
     cout << "\tdataType: DataExpanded" << endl;  
   }  
   if (dataType == -1) {  
     throw DataException("archiveData Error: undefined dataType");  
2232    }    }
2233      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2234    
2235    //      // Borrow DataTagged input from Data object
2236    // Collect data items common to all Data types      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2237    int noSamples = getNumSamples();      if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2238    int noDPPSample = getNumDataPointsPerSample();  
2239    int functionSpaceType = getFunctionSpace().getTypeCode();      // Prepare the DataConstant input
2240    int dataPointRank = getDataPointRank();      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2241    int dataPointSize = getDataPointSize();      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2242    int dataLength = getLength();  
2243    DataArrayView::ShapeType dataPointShape = getDataPointShape();      // Prepare a DataTagged output 2
2244    int referenceNumbers[noSamples];      res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2245    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      res.tag();
2246      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2247    }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2248    int tagNumbers[noSamples];  
2249    if (isTagged()) {      // Prepare offset into DataConstant
2250      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      int offset_1 = tmp_1->getPointOffset(0,0);
2251        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);      double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2252        // Get the views
2253    //     DataArrayView view_0 = tmp_0->getDefaultValue();
2254    //     DataArrayView view_2 = tmp_2->getDefaultValue();
2255    //     // Get the pointers to the actual data
2256    //     double *ptr_0 = &((view_0.getData())[0]);
2257    //     double *ptr_2 = &((view_2.getData())[0]);
2258    
2259        double *ptr_0 = &(tmp_0->getDefaultValue(0));
2260        double *ptr_2 = &(tmp_2->getDefaultValue(0));
2261    
2262        // Compute an MVP for the default
2263        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2264        // Compute an MVP for each tag
2265        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2266        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2267        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2268    //      tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2269    //       DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2270    //       DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2271    //       double *ptr_0 = &view_0.getData(0);
2272    //       double *ptr_2 = &view_2.getData(0);
2273    
2274          tmp_2->addTag(i->first);
2275          double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2276          double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2277          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2278        }
2279    
2280      }
2281      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2282    
2283        // Borrow DataTagged input from Data object
2284        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2285        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2286    
2287        // Borrow DataTagged input from Data object
2288        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2289        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2290    
2291        // Prepare a DataTagged output 2
2292        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2293        res.tag();  // DataTagged output
2294        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2295        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2296    
2297    //     // Get the views
2298    //     DataArrayView view_0 = tmp_0->getDefaultValue();
2299    //     DataArrayView view_1 = tmp_1->getDefaultValue();
2300    //     DataArrayView view_2 = tmp_2->getDefaultValue();
2301    //     // Get the pointers to the actual data
2302    //     double *ptr_0 = &((view_0.getData())[0]);
2303    //     double *ptr_1 = &((view_1.getData())[0]);
2304    //     double *ptr_2 = &((view_2.getData())[0]);
2305    
2306        double *ptr_0 = &(tmp_0->getDefaultValue(0));
2307        double *ptr_1 = &(tmp_1->getDefaultValue(0));
2308        double *ptr_2 = &(tmp_2->getDefaultValue(0));
2309    
2310    
2311        // Compute an MVP for the default
2312        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2313        // Merge the tags
2314        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2315        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2316        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2317        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2318          tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2319        }
2320        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2321          tmp_2->addTag(i->first);
2322        }
2323        // Compute an MVP for each tag
2324        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2325        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2326    //       DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2327    //       DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2328    //       DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2329    //       double *ptr_0 = &view_0.getData(0);
2330    //       double *ptr_1 = &view_1.getData(0);
2331    //       double *ptr_2 = &view_2.getData(0);
2332    
2333          double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2334          double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2335          double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2336    
2337          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2338        }
2339    
2340      }
2341      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2342    
2343        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2344        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2345        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2346        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2347        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2348        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2349        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2350        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2351        int sampleNo_0,dataPointNo_0;
2352        int numSamples_0 = arg_0_Z.getNumSamples();
2353        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2354        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2355        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2356          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2357          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2358          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2359            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2360            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2361            double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2362            double *ptr_2 = &(res.getDataAtOffset(offset_2));
2363            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2364          }
2365      }      }
   }  
   
   cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;  
   cout << "\tfunctionSpaceType: " << functionSpaceType << endl;  
   cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;  
2366    
   //  
   // Flatten Shape to an array of integers suitable for writing to file  
   int flatShape[4] = {0,0,0,0};  
   cout << "\tshape: < ";  
   for (int dim=0; dim<dataPointRank; dim++) {  
     flatShape[dim] = dataPointShape[dim];  
     cout << dataPointShape[dim] << " ";  
2367    }    }
2368    cout << ">" << endl;    else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2369    
2370    //      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2371    // Write common data items to archive file      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2372    ofstream archiveFile;      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2373    archiveFile.open(fileName.data(), ios::out);      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2374        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2375        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2376        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2377        int sampleNo_0,dataPointNo_0;
2378        int numSamples_0 = arg_0_Z.getNumSamples();
2379        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2380        int offset_1 = tmp_1->getPointOffset(0,0);
2381        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2382        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2383          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2384            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2385            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2386            double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2387            double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2388            double *ptr_2 = &(res.getDataAtOffset(offset_2));
2389            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2390          }
2391        }
2392    
   if (!archiveFile.good()) {  
     throw DataException("archiveData Error: problem opening archive file");  
   }  
2393    
   archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));  
   archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));  
   archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));  
   archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));  
   archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));  
   archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));  
   archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));  
   for (int dim = 0; dim < 4; dim++) {  
     archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));  
2394    }    }
2395    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2396      archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));  
2397    }      // After finding a common function space above the two inputs have the same numSamples and num DPPS
2398    if (isTagged()) {      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2399      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2400        archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));      DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2401        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2402        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2403        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2404        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2405        int sampleNo_0,dataPointNo_0;
2406        int numSamples_0 = arg_0_Z.getNumSamples();
2407        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2408        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2409        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2410          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2411          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2412          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2413            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2414            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2415            double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2416            double *ptr_2 = &(res.getDataAtOffset(offset_2));
2417            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2418          }
2419      }      }
   }  
2420    
   if (!archiveFile.good()) {  
     throw DataException("archiveData Error: problem writing to archive file");  
2421    }    }
2422      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2423    
2424    archiveFile.close();      // After finding a common function space above the two inputs have the same numSamples and num DPPS
2425        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2426        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2427        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2428        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2429        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2430        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2431        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2432        int sampleNo_0,dataPointNo_0;
2433        int numSamples_0 = arg_0_Z.getNumSamples();
2434        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2435        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2436        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2437          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2438            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2439            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2440            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2441            double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2442            double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2443            double *ptr_2 = &(res.getDataAtOffset(offset_2));
2444            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2445          }
2446        }
2447    
   if (!archiveFile.good()) {  
     throw DataException("archiveData Error: problem closing archive file");  
2448    }    }
2449      else {
2450    //      throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
   // Collect and archive underlying data values for each Data type  
   switch (dataType) {  
     case 0:  
       // DataEmpty  
       break;  
     case 1:  
       // DataConstant  
       break;  
     case 2:  
       // DataTagged  
       break;  
     case 3:  
       // DataExpanded  
       break;  
2451    }    }
2452    
2453      return res;
2454  }  }
2455    
2456  void  DataAbstract*
2457  Data::extractData(const std::string fileName,  Data::borrowData() const
                   const FunctionSpace& fspace)  
2458  {  {
2459    //    return m_data.get();
2460    // Can only extract Data to an object which is initially DataEmpty  }
   if (!isEmpty()) {  
     throw DataException("extractData Error: can only extract to DataEmpty object");  
   }  
2461    
   cout << "Extracting Data object from: " << fileName << endl;  
2462    
2463    int dataType;  std::string
2464    int noSamples;  Data::toString() const
2465    int noDPPSample;  {
2466    int functionSpaceType;      static const DataTypes::ValueType::size_type TOO_MANY_POINTS=80;
2467    int dataPointRank;      if (getNumDataPoints()*getDataPointSize()>TOO_MANY_POINTS)
2468    int dataPointSize;      {
2469    int dataLength;      stringstream temp;
2470    DataArrayView::ShapeType dataPointShape;      temp << "Summary: inf="<< inf() << " sup=" << sup() << " data points=" << getNumDataPoints();
2471    int flatShape[4];      return  temp.str();
2472        }
2473        return m_data->toString();
2474    }
2475    
   //  
   // Open the archive file and read common data items  
   ifstream archiveFile;  
   archiveFile.open(fileName.data(), ios::in);  
2476    
   if (!archiveFile.good()) {  
     throw DataException("extractData Error: problem opening archive file");  
   }  
2477    
2478    archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));  DataTypes::ValueType::const_reference
2479    archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));  Data::getDataAtOffset(DataTypes::ValueType::size_type i) const
2480    archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));  {
2481    archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));      return m_data->getDataAtOffset(i);
2482    archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));  }
   archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));  
   archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));  
   for (int dim = 0; dim < 4; dim++) {  
     archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));  
     if (flatShape[dim]>0) {  
       dataPointShape.push_back(flatShape[dim]);  
     }  
   }  
   int referenceNumbers[noSamples];  
   for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {  
     archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));  
   }  
   int tagNumbers[noSamples];  
   if (dataType==2) {  
     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {  
       archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));  
     }  
   }  
2483    
   if (!archiveFile.good()) {  
     throw DataException("extractData Error: problem reading from archive file");  
   }  
2484    
2485    archiveFile.close();  DataTypes::ValueType::reference
2486    Data::getDataAtOffset(DataTypes::ValueType::size_type i)
2487    {
2488        return m_data->getDataAtOffset(i);
2489    }
2490    
2491    if (!archiveFile.good()) {  DataTypes::ValueType::const_reference
2492      throw DataException("extractData Error: problem closing archive file");  Data::getDataPoint(int sampleNo, int dataPointNo) const
2493    }  {
2494        return m_data->getDataAtOffset(m_data->getPointOffset(sampleNo, dataPointNo));
2495    }
2496    
   switch (dataType) {  
     case 0:  
       cout << "\tdataType: DataEmpty" << endl;  
       break;  
     case 1:  
       cout << "\tdataType: DataConstant" << endl;  
       break;  
     case 2:  
       cout << "\tdataType: DataTagged" << endl;  
       break;  
     case 3:  
       cout << "\tdataType: DataExpanded" << endl;  
       break;  
     default:  
       throw DataException("extractData Error: undefined dataType read from archive file");  
       break;  
   }  
2497    
2498    cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;  DataTypes::ValueType::reference
2499    cout << "\tfunctionSpaceType: " << functionSpaceType << endl;  Data::getDataPoint(int sampleNo, int dataPointNo)
2500    cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;  {
2501    cout << "\tshape: < ";      return m_data->getDataAtOffset(m_data->getPointOffset(sampleNo, dataPointNo));
2502    for (int dim = 0; dim < dataPointRank; dim++) {  }
     cout << dataPointShape[dim] << " ";  
   }  
   cout << ">" << endl;  
2503    
   //  
   // Verify that supplied FunctionSpace object is compatible with this Data object.  
   if ( (fspace.getTypeCode()!=functionSpaceType) ||  
        (fspace.getNumSamples()!=noSamples) ||  
        (fspace.getNumDPPSample()!=noDPPSample)  
      ) {  
     throw DataException("extractData Error: incompatible FunctionSpace");  
   }  
   for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {  
     if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {  
       throw DataException("extractData Error: incompatible FunctionSpace");  
     }  
   }  
   if (dataType==2) {  
     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {  
       if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {  
         throw DataException("extractData Error: incompatible FunctionSpace");  
       }  
     }  
   }  
2504    
2505    //  /* Member functions specific to the MPI implementation */
   // Construct a DataVector to hold underlying data values  
   DataVector dataVec(dataLength);  
2506    
2507    //  void
2508    // Load this DataVector with the appropriate values  Data::print()
2509    switch (dataType) {  {
2510      case 0:    int i,j;
       // DataEmpty  
       break;  
     case 1:  
       // DataConstant  
       break;  
     case 2:  
       // DataTagged  
       break;  
     case 3:  
       // DataExpanded  
       break;  
   }  
2511    
2512    //    printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2513    // Construct an appropriate Data object    for( i=0; i<getNumSamples(); i++ )
2514    DataAbstract* tempData;    {
2515    switch (dataType) {      printf( "[%6d]", i );
2516      case 0:      for( j=0; j<getNumDataPointsPerSample(); j++ )
2517        // DataEmpty        printf( "\t%10.7g", (getSampleData(i))[j] );
2518        tempData=new DataEmpty();      printf( "\n" );
       break;  
     case 1:  
       // DataConstant  
       tempData=new DataConstant(fspace,dataPointShape,dataVec);  
       break;  
     case 2:  
       // DataTagged  
       tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);  
       break;  
     case 3:  
       // DataExpanded  
       tempData=new DataExpanded(fspace,dataPointShape,dataVec);  
       break;  
2519    }    }
2520    shared_ptr<DataAbstract> temp_data(tempData);  }
2521    m_data=temp_data;  void
2522    Data::dump(const std::string fileName) const
2523    {
2524      try
2525         {
2526            return m_data->dump(fileName);
2527         }
2528         catch (exception& e)
2529         {
2530            cout << e.what() << endl;
2531         }
2532    }
2533    
2534    int
2535    Data::get_MPISize() const
2536    {
2537        int size;
2538    #ifdef PASO_MPI
2539        int error;
2540        error = MPI_Comm_size( get_MPIComm(), &size );
2541    #else
2542        size = 1;
2543    #endif
2544        return size;
2545  }  }
2546    
2547  ostream& escript::operator<<(ostream& o, const Data& data)  int
2548    Data::get_MPIRank() const
2549  {  {
2550    o << data.toString();      int rank;
2551    return o;  #ifdef PASO_MPI
2552        int error;
2553        error = MPI_Comm_rank( get_MPIComm(), &rank );
2554    #else
2555        rank = 0;
2556    #endif
2557        return rank;
2558    }
2559    
2560    MPI_Comm
2561    Data::get_MPIComm() const
2562    {
2563    #ifdef PASO_MPI
2564        return MPI_COMM_WORLD;
2565    #else
2566        return -1;
2567    #endif
2568  }  }
2569    
2570    

Legend:
Removed from v.121  
changed lines
  Added in v.1872

  ViewVC Help
Powered by ViewVC 1.1.26