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

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

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

revision 2037 by jfenwick, Thu Nov 13 06:17:12 2008 UTC revision 2673 by jfenwick, Fri Sep 18 05:33:10 2009 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2009 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 26  Line 26 
26  #include "EscriptParams.h"  #include "EscriptParams.h"
27    
28  extern "C" {  extern "C" {
29  #include "escript/blocktimer.h"  #include "esysUtils/blocktimer.h"
30  }  }
31    
32  #include <fstream>  #include <fstream>
33  #include <algorithm>  #include <algorithm>
34  #include <vector>  #include <vector>
35  #include <functional>  #include <functional>
36    #include <sstream>  // so we can throw messages about ranks
37    
38  #include <boost/python/dict.hpp>  #include <boost/python/dict.hpp>
39  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
40  #include <boost/python/long.hpp>  #include <boost/python/long.hpp>
41    #include "WrappedArray.h"
42    
43  using namespace std;  using namespace std;
44  using namespace boost::python;  using namespace boost::python;
# Line 45  using namespace escript; Line 47  using namespace escript;
47    
48  // ensure the current object is not a DataLazy  // ensure the current object is not a DataLazy
49  // The idea was that we could add an optional warning whenever a resolve is forced  // The idea was that we could add an optional warning whenever a resolve is forced
50  #define FORCERESOLVE if (isLazy()) {resolve();}  // #define forceResolve() if (isLazy()) {#resolve();}
51    
52    #define AUTOLAZYON escriptParams.getAUTOLAZY()
53    #define MAKELAZYOP(X)   if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
54      {\
55        DataLazy* c=new DataLazy(borrowDataPtr(),X);\
56        return Data(c);\
57      }
58    #define MAKELAZYOPOFF(X,Y) if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
59      {\
60        DataLazy* c=new DataLazy(borrowDataPtr(),X,Y);\
61        return Data(c);\
62      }
63    
64    #define MAKELAZYOP2(X,Y,Z) if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
65      {\
66        DataLazy* c=new DataLazy(borrowDataPtr(),X,Y,Z);\
67        return Data(c);\
68      }
69    
70    #define MAKELAZYBINSELF(R,X)   if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
71      {\
72        DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
73    /*         m_data=c->getPtr();*/     set_m_data(c->getPtr());\
74        return (*this);\
75      }
76    
77    // like the above but returns a new data rather than *this
78    #define MAKELAZYBIN(R,X)   if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
79      {\
80        DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
81        return Data(c);\
82      }
83    
84    #define MAKELAZYBIN2(L,R,X)   if (L.isLazy() || R.isLazy() || (AUTOLAZYON && (L.isExpanded() || R.isExpanded()))) \
85      {\
86        DataLazy* c=new DataLazy(L.borrowDataPtr(),R.borrowDataPtr(),X);\
87        return Data(c);\
88      }
89    
90    namespace
91    {
92    
93    template <class ARR>
94    inline
95    boost::python::tuple
96    pointToTuple1(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
97    {
98        using namespace boost::python;
99        using boost::python::tuple;
100        using boost::python::list;
101    
102        list l;
103        unsigned int dim=shape[0];
104        for (size_t i=0;i<dim;++i)
105        {
106        l.append(v[i+offset]);
107        }
108        return tuple(l);
109    }
110    
111    template <class ARR>
112    inline
113    boost::python::tuple
114    pointToTuple2(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
115    {
116        using namespace boost::python;
117        using boost::python::tuple;
118        using boost::python::list;
119    
120        unsigned int shape0=shape[0];
121        unsigned int shape1=shape[1];
122        list lj;
123        for (size_t j=0;j<shape0;++j)
124        {
125            list li;
126        for (size_t i=0;i<shape1;++i)
127        {
128            li.append(v[offset+DataTypes::getRelIndex(shape,j,i)]);
129        }
130        lj.append(tuple(li));
131        }
132        return tuple(lj);
133    }
134    
135    template <class ARR>
136    inline
137    boost::python::tuple
138    pointToTuple3(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
139    {
140        using namespace boost::python;
141        using boost::python::tuple;
142        using boost::python::list;
143    
144        unsigned int shape0=shape[0];
145        unsigned int shape1=shape[1];
146        unsigned int shape2=shape[2];
147    
148        list lk;
149        for (size_t k=0;k<shape0;++k)
150        {
151            list lj;
152        for (size_t j=0;j<shape1;++j)
153        {
154            list li;
155            for (size_t i=0;i<shape2;++i)
156            {
157                    li.append(v[offset+DataTypes::getRelIndex(shape,k,j,i)]);
158                }
159            lj.append(tuple(li));
160            }
161            lk.append(tuple(lj));
162        }
163        return tuple(lk);
164    }
165    
166    template <class ARR>
167    inline
168    boost::python::tuple
169    pointToTuple4(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
170    {
171        using namespace boost::python;
172        using boost::python::tuple;
173        using boost::python::list;
174    
175        unsigned int shape0=shape[0];
176        unsigned int shape1=shape[1];
177        unsigned int shape2=shape[2];
178        unsigned int shape3=shape[3];
179    
180        list ll;
181        for (size_t l=0;l<shape0;++l)
182        {
183            list lk;
184        for (size_t k=0;k<shape1;++k)
185        {
186                list lj;
187                for (size_t j=0;j<shape2;++j)
188                {
189                    list li;
190                    for (size_t i=0;i<shape3;++i)
191                    {
192                        li.append(v[offset+DataTypes::getRelIndex(shape,l,k,j,i)]);
193                    }
194                    lj.append(tuple(li));
195                }
196                lk.append(tuple(lj));
197        }
198            ll.append(tuple(lk));
199        }
200        return tuple(ll);
201    }
202    
203    
204    // This should be safer once the DataC RO changes have been brought in
205    template <class ARR>
206    boost::python::tuple
207    pointToTuple( const DataTypes::ShapeType& shape,ARR v)
208    {
209       using namespace boost::python;
210       using boost::python::list;
211       int rank=shape.size();
212       if (rank==0)
213       {
214        return make_tuple(v[0]);
215       }
216       else if (rank==1)
217       {
218            return pointToTuple1(shape,v,0);
219       }
220       else if (rank==2)
221       {
222        return pointToTuple2(shape,v,0);
223       }
224       else if (rank==3)
225       {
226        return pointToTuple3(shape,v,0);
227       }
228       else if (rank==4)
229       {
230        return pointToTuple4(shape,v,0);
231       }
232       else
233         throw DataException("Unknown rank in pointToTuple.");
234    }
235    
236    }  // anonymous namespace
237    
238  Data::Data()  Data::Data()
239        : m_shared(false), m_lazy(false)
240  {  {
241    //    //
242    // Default data is type DataEmpty    // Default data is type DataEmpty
243    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
244    m_data=temp->getPtr();  //   m_data=temp->getPtr();
245      set_m_data(temp->getPtr());
246    m_protected=false;    m_protected=false;
247  }  }
248    
# Line 60  Data::Data(double value, Line 250  Data::Data(double value,
250             const tuple& shape,             const tuple& shape,
251             const FunctionSpace& what,             const FunctionSpace& what,
252             bool expanded)             bool expanded)
253        : m_shared(false), m_lazy(false)
254  {  {
255    DataTypes::ShapeType dataPointShape;    DataTypes::ShapeType dataPointShape;
256    for (int i = 0; i < shape.attr("__len__")(); ++i) {    for (int i = 0; i < shape.attr("__len__")(); ++i) {
# Line 76  Data::Data(double value, Line 267  Data::Data(double value,
267         const DataTypes::ShapeType& dataPointShape,         const DataTypes::ShapeType& dataPointShape,
268         const FunctionSpace& what,         const FunctionSpace& what,
269             bool expanded)             bool expanded)
270        : m_shared(false), m_lazy(false)
271  {  {
272    int len = DataTypes::noValues(dataPointShape);    int len = DataTypes::noValues(dataPointShape);
   
273    DataVector temp_data(len,value,len);    DataVector temp_data(len,value,len);
 //   DataArrayView temp_dataView(temp_data, dataPointShape);  
   
 //   initialise(temp_dataView, what, expanded);  
274    initialise(temp_data, dataPointShape, what, expanded);    initialise(temp_data, dataPointShape, what, expanded);
   
275    m_protected=false;    m_protected=false;
276  }  }
277    
278  Data::Data(const Data& inData)  Data::Data(const Data& inData)
279        : m_shared(false), m_lazy(false)
280  {  {
281    m_data=inData.m_data;    set_m_data(inData.m_data);
282    m_protected=inData.isProtected();    m_protected=inData.isProtected();
283  }  }
284    
285    
286  Data::Data(const Data& inData,  Data::Data(const Data& inData,
287             const DataTypes::RegionType& region)             const DataTypes::RegionType& region)
288        : m_shared(false), m_lazy(false)
289  {  {
290    DataAbstract_ptr dat=inData.m_data;    DataAbstract_ptr dat=inData.m_data;
291    if (inData.isLazy())    if (inData.isLazy())
# Line 110  Data::Data(const Data& inData, Line 299  Data::Data(const Data& inData,
299    //    //
300    // Create Data which is a slice of another Data    // Create Data which is a slice of another Data
301    DataAbstract* tmp = dat->getSlice(region);    DataAbstract* tmp = dat->getSlice(region);
302    m_data=DataAbstract_ptr(tmp);    set_m_data(DataAbstract_ptr(tmp));
303    m_protected=false;    m_protected=false;
304    
305  }  }
306    
307  Data::Data(const Data& inData,  Data::Data(const Data& inData,
308             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
309        : m_shared(false), m_lazy(false)
310  {  {
311    if (inData.isEmpty())    if (inData.isEmpty())
312    {    {
313      throw DataException("Error - will not interpolate for instances of DataEmpty.");      throw DataException("Error - will not interpolate for instances of DataEmpty.");
314    }    }
315    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
316      m_data=inData.m_data;      set_m_data(inData.m_data);
317    }    }
318    else    else
319    {    {
# Line 130  Data::Data(const Data& inData, Line 321  Data::Data(const Data& inData,
321      if (inData.isConstant()) {  // for a constant function, we just need to use the new function space      if (inData.isConstant()) {  // for a constant function, we just need to use the new function space
322        if (!inData.probeInterpolation(functionspace))        if (!inData.probeInterpolation(functionspace))
323        {           // Even though this is constant, we still need to check whether interpolation is allowed        {           // Even though this is constant, we still need to check whether interpolation is allowed
324      throw FunctionSpaceException("Call to probeInterpolation returned false for DataConstant.");      throw FunctionSpaceException("Cannot interpolate across to the domain of the specified FunctionSpace. (DataConstant)");
325        }        }
326        // if the data is not lazy, this will just be a cast to DataReady        // if the data is not lazy, this will just be a cast to DataReady
327        DataReady_ptr dr=inData.m_data->resolve();        DataReady_ptr dr=inData.m_data->resolve();
328        DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVector());          DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVectorRO());
329        m_data=DataAbstract_ptr(dc);  //       m_data=DataAbstract_ptr(dc);
330          set_m_data(DataAbstract_ptr(dc));
331      } else {      } else {
332        Data tmp(0,inData.getDataPointShape(),functionspace,true);        Data tmp(0,inData.getDataPointShape(),functionspace,true);
333        // Note: Must use a reference or pointer to a derived object        // Note: Must use a reference or pointer to a derived object
# Line 149  Data::Data(const Data& inData, Line 341  Data::Data(const Data& inData,
341        } else {        } else {
342          inDataDomain->interpolateACross(tmp,inData);          inDataDomain->interpolateACross(tmp,inData);
343        }        }
344        m_data=tmp.m_data;  //       m_data=tmp.m_data;
345          set_m_data(tmp.m_data);
346      }      }
347    }    }
348    m_protected=false;    m_protected=false;
349  }  }
350    
351  Data::Data(DataAbstract* underlyingdata)  Data::Data(DataAbstract* underlyingdata)
352        : m_shared(false), m_lazy(false)
353  {  {
354  //  m_data=shared_ptr<DataAbstract>(underlyingdata);      set_m_data(underlyingdata->getPtr());
     m_data=underlyingdata->getPtr();  
355      m_protected=false;      m_protected=false;
356  }  }
357    
358  Data::Data(DataAbstract_ptr underlyingdata)  Data::Data(DataAbstract_ptr underlyingdata)
359        : m_shared(false), m_lazy(false)
360  {  {
361      m_data=underlyingdata;      set_m_data(underlyingdata);
362      m_protected=false;      m_protected=false;
363  }  }
364    
   
 Data::Data(const numeric::array& value,  
        const FunctionSpace& what,  
            bool expanded)  
 {  
   initialise(value,what,expanded);  
   m_protected=false;  
 }  
 /*  
 Data::Data(const DataArrayView& value,  
        const FunctionSpace& what,  
            bool expanded)  
 {  
   initialise(value,what,expanded);  
   m_protected=false;  
 }*/  
   
365  Data::Data(const DataTypes::ValueType& value,  Data::Data(const DataTypes::ValueType& value,
366           const DataTypes::ShapeType& shape,           const DataTypes::ShapeType& shape,
367                   const FunctionSpace& what,                   const FunctionSpace& what,
368                   bool expanded)                   bool expanded)
369        : m_shared(false), m_lazy(false)
370  {  {
371     initialise(value,shape,what,expanded);     initialise(value,shape,what,expanded);
372     m_protected=false;     m_protected=false;
# Line 198  Data::Data(const DataTypes::ValueType& v Line 376  Data::Data(const DataTypes::ValueType& v
376  Data::Data(const object& value,  Data::Data(const object& value,
377         const FunctionSpace& what,         const FunctionSpace& what,
378             bool expanded)             bool expanded)
379        : m_shared(false), m_lazy(false)
380  {  {
381    numeric::array asNumArray(value);    WrappedArray w(value);
382    initialise(asNumArray,what,expanded);    initialise(w,what,expanded);
383    m_protected=false;    m_protected=false;
384  }  }
385    
386    
387  Data::Data(const object& value,  Data::Data(const object& value,
388             const Data& other)             const Data& other)
389        : m_shared(false), m_lazy(false)
390  {  {
391    numeric::array asNumArray(value);    WrappedArray w(value);
392    
393    // extract the shape of the numarray    // extract the shape of the array
394    DataTypes::ShapeType tempShape=DataTypes::shapeFromNumArray(asNumArray);    const DataTypes::ShapeType& tempShape=w.getShape();
395  // /*  for (int i=0; i < asNumArray.getrank(); i++) {    if (w.getRank()==0) {
 //     tempShape.push_back(extract<int>(asNumArray.getshape()[i]));  
 //   }*/  
 //   // get the space for the data vector  
 //   int len = DataTypes::noValues(tempShape);  
 //   DataVector temp_data(len, 0.0, len);  
 // /*  DataArrayView temp_dataView(temp_data, tempShape);  
 //   temp_dataView.copy(asNumArray);*/  
 //   temp_data.copyFromNumArray(asNumArray);  
   
   //  
   // Create DataConstant using the given value and all other parameters  
   // copied from other. If value is a rank 0 object this Data  
   // will assume the point data shape of other.  
   
   if (DataTypes::getRank(tempShape)/*temp_dataView.getRank()*/==0) {  
396    
397    
398      // get the space for the data vector      // get the space for the data vector
399      int len1 = DataTypes::noValues(tempShape);      int len1 = DataTypes::noValues(tempShape);
400      DataVector temp_data(len1, 0.0, len1);      DataVector temp_data(len1, 0.0, len1);
401      temp_data.copyFromNumArray(asNumArray);      temp_data.copyFromArray(w,1);
402    
403      int len = DataTypes::noValues(other.getDataPointShape());      int len = DataTypes::noValues(other.getDataPointShape());
404    
405      DataVector temp2_data(len, temp_data[0]/*temp_dataView()*/, len);      DataVector temp2_data(len, temp_data[0], len);
     //DataArrayView temp2_dataView(temp2_data, other.getPointDataView().getShape());  
 //     initialise(temp2_dataView, other.getFunctionSpace(), false);  
   
406      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
407  //     boost::shared_ptr<DataAbstract> sp(t);  //     m_data=DataAbstract_ptr(t);
408  //     m_data=sp;      set_m_data(DataAbstract_ptr(t));
     m_data=DataAbstract_ptr(t);  
409    
410    } else {    } else {
411      //      //
412      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
413  //     initialise(temp_dataView, other.getFunctionSpace(), false);      DataConstant* t=new DataConstant(w,other.getFunctionSpace());
414      DataConstant* t=new DataConstant(asNumArray,other.getFunctionSpace());  //     m_data=DataAbstract_ptr(t);
415  //     boost::shared_ptr<DataAbstract> sp(t);      set_m_data(DataAbstract_ptr(t));
 //     m_data=sp;  
     m_data=DataAbstract_ptr(t);  
416    }    }
417    m_protected=false;    m_protected=false;
418  }  }
419    
420  Data::~Data()  Data::~Data()
421  {  {
422      set_m_data(DataAbstract_ptr());
423  }  }
424    
425    
426    // only call in thread safe contexts.
427    // This method should be atomic
428    void Data::set_m_data(DataAbstract_ptr p)
429    {
430      if (m_data.get()!=0)  // release old ownership
431      {
432        m_data->removeOwner(this);
433      }
434      if (p.get()!=0)
435      {
436        m_data=p;
437        m_data->addOwner(this);
438        m_shared=m_data->isShared();
439        m_lazy=m_data->isLazy();
440      }
441    }
442    
443  void  void Data::initialise(const WrappedArray& value,
 Data::initialise(const boost::python::numeric::array& value,  
444                   const FunctionSpace& what,                   const FunctionSpace& what,
445                   bool expanded)                   bool expanded)
446  {  {
# Line 277  Data::initialise(const boost::python::nu Line 451  Data::initialise(const boost::python::nu
451    // within the shared_ptr constructor.    // within the shared_ptr constructor.
452    if (expanded) {    if (expanded) {
453      DataAbstract* temp=new DataExpanded(value, what);      DataAbstract* temp=new DataExpanded(value, what);
454  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
455  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
456    } else {    } else {
457      DataAbstract* temp=new DataConstant(value, what);      DataAbstract* temp=new DataConstant(value, what);
458  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
459  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
460    }    }
461  }  }
462    
# Line 302  Data::initialise(const DataTypes::ValueT Line 474  Data::initialise(const DataTypes::ValueT
474    // within the shared_ptr constructor.    // within the shared_ptr constructor.
475    if (expanded) {    if (expanded) {
476      DataAbstract* temp=new DataExpanded(what, shape, value);      DataAbstract* temp=new DataExpanded(what, shape, value);
477  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
478  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
479    } else {    } else {
480      DataAbstract* temp=new DataConstant(what, shape, value);      DataAbstract* temp=new DataConstant(what, shape, value);
481  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
482  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
483    }    }
484  }  }
485    
486    
 // void  
 // Data::CompareDebug(const Data& rd)  
 // {  
 //  using namespace std;  
 //  bool mismatch=false;  
 //  std::cout << "Comparing left and right" << endl;  
 //  const DataTagged* left=dynamic_cast<DataTagged*>(m_data.get());  
 //  const DataTagged* right=dynamic_cast<DataTagged*>(rd.m_data.get());  
 //    
 //  if (left==0)  
 //  {  
 //      cout << "left arg is not a DataTagged\n";  
 //      return;  
 //  }  
 //    
 //  if (right==0)  
 //  {  
 //      cout << "right arg is not a DataTagged\n";  
 //      return;  
 //  }  
 //  cout << "Num elements=" << left->getVector().size() << ":" << right->getVector().size() << std::endl;  
 //  cout << "Shapes ";  
 //  if (left->getShape()==right->getShape())  
 //  {  
 //      cout << "ok\n";  
 //  }  
 //  else  
 //  {  
 //      cout << "Problem: shapes do not match\n";  
 //      mismatch=true;  
 //  }  
 //  int lim=left->getVector().size();  
 //  if (right->getVector().size()) lim=right->getVector().size();  
 //  for (int i=0;i<lim;++i)  
 //  {  
 //      if (left->getVector()[i]!=right->getVector()[i])  
 //      {  
 //          cout << "[" << i << "] value mismatch " << left->getVector()[i] << ":" << right->getVector()[i] << endl;  
 //          mismatch=true;  
 //      }  
 //  }  
 //  
 //  // still need to check the tag map  
 //  // also need to watch what is happening to function spaces, are they copied or what?  
 //  
 //  const DataTagged::DataMapType& mapleft=left->getTagLookup();  
 //  const DataTagged::DataMapType& mapright=right->getTagLookup();  
 //  
 //  if (mapleft.size()!=mapright.size())  
 //  {  
 //      cout << "Maps are different sizes " << mapleft.size() << ":" << mapright.size() << endl;  
 //      mismatch=true;  
 //      cout << "Left map\n";  
 //      DataTagged::DataMapType::const_iterator i,j;  
 //      for (i=mapleft.begin();i!=mapleft.end();++i) {  
 //          cout << "(" << i->first << "=>" << i->second << ")\n";  
 //      }  
 //      cout << "Right map\n";  
 //      for (i=mapright.begin();i!=mapright.end();++i) {  
 //          cout << "(" << i->first << "=>" << i->second << ")\n";  
 //      }  
 //      cout << "End map\n";  
 //  
 //  }  
 //  
 //  DataTagged::DataMapType::const_iterator i,j;  
 //  for (i=mapleft.begin(),j=mapright.begin();i!=mapleft.end() && j!=mapright.end();++i,++j) {  
 //     if ((i->first!=j->first) || (i->second!=j->second))  
 //     {  
 //      cout << "(" << i->first << "=>" << i->second << ")";  
 //      cout << ":(" << j->first << "=>" << j->second << ") ";  
 //      mismatch=true;  
 //            }  
 //  }  
 //  if (mismatch)  
 //  {  
 //      cout << "#Mismatch\n";  
 //  }  
 // }  
   
487  escriptDataC  escriptDataC
488  Data::getDataC()  Data::getDataC()
489  {  {
# Line 410  Data::getDataC() const Line 500  Data::getDataC() const
500    return temp;    return temp;
501  }  }
502    
503    size_t
504    Data::getSampleBufferSize() const
505    {
506      return m_data->getSampleBufferSize();
507    }
508    
509    
510  const boost::python::tuple  const boost::python::tuple
511  Data::getShapeTuple() const  Data::getShapeTuple() const
512  {  {
# Line 435  Data::getShapeTuple() const Line 532  Data::getShapeTuple() const
532  // It can't work out what type the function is based soley on its name.  // It can't work out what type the function is based soley on its name.
533  // There are ways to fix this involving creating function pointer variables for each form  // There are ways to fix this involving creating function pointer variables for each form
534  // but there doesn't seem to be a need given that the methods have the same name from the python point of view  // but there doesn't seem to be a need given that the methods have the same name from the python point of view
535  Data*  Data
536  Data::copySelf()  Data::copySelf()
537  {  {
538     DataAbstract* temp=m_data->deepCopy();     DataAbstract* temp=m_data->deepCopy();
539     return new Data(temp);     return Data(temp);
540  }  }
541    
542  void  void
# Line 447  Data::copy(const Data& other) Line 544  Data::copy(const Data& other)
544  {  {
545    DataAbstract* temp=other.m_data->deepCopy();    DataAbstract* temp=other.m_data->deepCopy();
546    DataAbstract_ptr p=temp->getPtr();    DataAbstract_ptr p=temp->getPtr();
547    m_data=p;  //   m_data=p;
548      set_m_data(p);
549  }  }
550    
551    
# Line 463  Data::delaySelf() Line 561  Data::delaySelf()
561  {  {
562    if (!isLazy())    if (!isLazy())
563    {    {
564      m_data=(new DataLazy(m_data))->getPtr();  //  m_data=(new DataLazy(m_data))->getPtr();
565        set_m_data((new DataLazy(m_data))->getPtr());
566    }    }
567  }  }
568    
569    
570    // For lazy data, it would seem that DataTagged will need to be treated differently since even after setting all tags
571    // to zero, all the tags from all the DataTags would be in the result.
572    // However since they all have the same value (0) whether they are there or not should not matter.
573    // So I have decided that for all types this method will create a constant 0.
574    // It can be promoted up as required.
575    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
576    // but we can deal with that if it arises.
577    //
578  void  void
579  Data::setToZero()  Data::setToZero()
580  {  {
# Line 474  Data::setToZero() Line 582  Data::setToZero()
582    {    {
583       throw DataException("Error - Operations not permitted on instances of DataEmpty.");       throw DataException("Error - Operations not permitted on instances of DataEmpty.");
584    }    }
585    m_data->setToZero();    if (isLazy())
586      {
587         DataTypes::ValueType v(getNoValues(),0);
588         DataConstant* dc=new DataConstant(getFunctionSpace(),getDataPointShape(),v);
589         DataLazy* dl=new DataLazy(dc->getPtr());
590         set_m_data(dl->getPtr());
591      }
592      else
593      {
594         exclusiveWrite();
595         m_data->setToZero();
596      }
597  }  }
598    
599    
600  void  void
601  Data::copyWithMask(const Data& other,  Data::copyWithMask(const Data& other,
602                     const Data& mask)                     const Data& mask)
# Line 538  Data::copyWithMask(const Data& other, Line 658  Data::copyWithMask(const Data& other,
658    {    {
659      throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");      throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
660    }    }
661      unsigned int selfrank=getDataPointRank();
662      unsigned int otherrank=other2.getDataPointRank();
663      unsigned int maskrank=mask2.getDataPointRank();
664      if ((selfrank==0) && (otherrank>0 || maskrank>0))
665      {
666        // to get here we must be copying from a large object into a scalar
667        // I am not allowing this.
668        // If you are calling copyWithMask then you are considering keeping some existing values
669        // and so I'm going to assume that you don't want your data objects getting a new shape.
670        throw DataException("Attempt to copyWithMask into a scalar from an object or mask with rank>0.");
671      }
672      exclusiveWrite();
673    // Now we iterate over the elements    // Now we iterate over the elements
674    DataVector& self=getReadyPtr()->getVector();    DataVector& self=getReady()->getVectorRW();;
675    const DataVector& ovec=other2.getReadyPtr()->getVector();    const DataVector& ovec=other2.getReadyPtr()->getVectorRO();
676    const DataVector& mvec=mask2.getReadyPtr()->getVector();    const DataVector& mvec=mask2.getReadyPtr()->getVectorRO();
677    if ((self.size()!=ovec.size()) || (self.size()!=mvec.size()))  
678      if ((selfrank>0) && (otherrank==0) &&(maskrank==0))
679      {
680        // Not allowing this combination.
681        // it is not clear what the rank of the target should be.
682        // Should it be filled with the scalar (rank stays the same);
683        // or should the target object be reshaped to be a scalar as well.
684        throw DataException("Attempt to copyWithMask from scalar mask and data into non-scalar target.");
685      }
686      if ((selfrank>0) && (otherrank>0) &&(maskrank==0))
687    {    {
688      throw DataException("Error - size mismatch in arguments to copyWithMask.");      if (mvec[0]>0)      // copy whole object if scalar is >0
689        {
690            copy(other);
691        }
692        return;
693      }
694      if (isTagged())       // so all objects involved will also be tagged
695      {
696        // note the !
697        if (!((getDataPointShape()==mask2.getDataPointShape()) &&
698            ((other2.getDataPointShape()==mask2.getDataPointShape()) || (otherrank==0))))
699        {
700            throw DataException("copyWithMask, shape mismatch.");
701        }
702    
703        // We need to consider the possibility that tags are missing or in the wrong order
704        // My guiding assumption here is: All tagged Datas are assumed to have the default value for
705        // all tags which are not explicitly defined
706    
707        const DataTagged* mptr=dynamic_cast<const DataTagged*>(mask2.m_data.get());
708        const DataTagged* optr=dynamic_cast<const DataTagged*>(other2.m_data.get());
709        DataTagged* tptr=dynamic_cast<DataTagged*>(m_data.get());
710    
711        // first, add any tags missing from other or mask
712        const DataTagged::DataMapType& olookup=optr->getTagLookup();
713            const DataTagged::DataMapType& mlookup=mptr->getTagLookup();
714        const DataTagged::DataMapType& tlookup=tptr->getTagLookup();
715        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
716        for (i=olookup.begin();i!=olookup.end();i++)
717        {
718               tptr->addTag(i->first);
719            }
720            for (i=mlookup.begin();i!=mlookup.end();i++) {
721               tptr->addTag(i->first);
722            }
723        // now we know that *this has all the required tags but they aren't guaranteed to be in
724        // the same order
725    
726        // There are two possibilities: 1. all objects have the same rank. 2. other is a scalar
727        if ((selfrank==otherrank) && (otherrank==maskrank))
728        {
729            for (i=tlookup.begin();i!=tlookup.end();i++)
730            {
731                // get the target offset
732                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
733                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
734                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
735                for (int j=0;j<getDataPointSize();++j)
736                {
737                    if (mvec[j+moff]>0)
738                    {
739                        self[j+toff]=ovec[j+ooff];
740                    }
741                }
742                }
743            // now for the default value
744            for (int j=0;j<getDataPointSize();++j)
745            {
746                if (mvec[j+mptr->getDefaultOffset()]>0)
747                {
748                    self[j+tptr->getDefaultOffset()]=ovec[j+optr->getDefaultOffset()];
749                }
750            }
751        }
752        else    // other is a scalar
753        {
754            for (i=tlookup.begin();i!=tlookup.end();i++)
755            {
756                // get the target offset
757                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
758                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
759                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
760                for (int j=0;j<getDataPointSize();++j)
761                {
762                    if (mvec[j+moff]>0)
763                    {
764                        self[j+toff]=ovec[ooff];
765                    }
766                }
767                }
768            // now for the default value
769            for (int j=0;j<getDataPointSize();++j)
770            {
771                if (mvec[j+mptr->getDefaultOffset()]>0)
772                {
773                    self[j+tptr->getDefaultOffset()]=ovec[0];
774                }
775            }
776        }
777    
778        return;         // ugly
779      }
780      // mixed scalar and non-scalar operation
781      if ((selfrank>0) && (otherrank==0) && (mask2.getDataPointShape()==getDataPointShape()))
782      {
783            size_t num_points=self.size();
784        // OPENMP 3.0 allows unsigned loop vars.
785        #if defined(_OPENMP) && (_OPENMP < 200805)
786        long i;
787        #else
788        size_t i;
789        #endif  
790        size_t psize=getDataPointSize();    
791        #pragma omp parallel for private(i) schedule(static)
792        for (i=0;i<num_points;++i)
793        {
794            if (mvec[i]>0)
795            {
796                self[i]=ovec[i/psize];      // since this is expanded there is one scalar
797            }                   // dest point
798        }
799        return;         // ugly!
800      }
801      // tagged data is already taken care of so we only need to worry about shapes
802      // special cases with scalars are already dealt with so all we need to worry about is shape
803      if ((getDataPointShape()!=other2.getDataPointShape()) || getDataPointShape()!=mask2.getDataPointShape())
804      {
805        ostringstream oss;
806        oss <<"Error - size mismatch in arguments to copyWithMask.";
807        oss << "\nself_shape=" << DataTypes::shapeToString(getDataPointShape());
808        oss << " other2_shape=" << DataTypes::shapeToString(other2.getDataPointShape());
809        oss << " mask2_shape=" << DataTypes::shapeToString(mask2.getDataPointShape());
810        throw DataException(oss.str());
811    }    }
812    size_t num_points=self.size();    size_t num_points=self.size();
813    
# Line 564  Data::copyWithMask(const Data& other, Line 827  Data::copyWithMask(const Data& other,
827    }    }
828  }  }
829    
   
   
830  bool  bool
831  Data::isExpanded() const  Data::isExpanded() const
832  {  {
# Line 574  Data::isExpanded() const Line 835  Data::isExpanded() const
835  }  }
836    
837  bool  bool
838    Data::actsExpanded() const
839    {
840      return m_data->actsExpanded();
841    
842    }
843    
844    
845    bool
846  Data::isTagged() const  Data::isTagged() const
847  {  {
848    DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());    DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
# Line 597  Data::isConstant() const Line 866  Data::isConstant() const
866  bool  bool
867  Data::isLazy() const  Data::isLazy() const
868  {  {
869    return m_data->isLazy();    return m_lazy;    // not asking m_data because we need to be able to ask this while m_data is changing
870  }  }
871    
872  // at the moment this is synonymous with !isLazy() but that could change  // at the moment this is synonymous with !isLazy() but that could change
# Line 628  Data::expand() Line 897  Data::expand()
897    if (isConstant()) {    if (isConstant()) {
898      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
899      DataAbstract* temp=new DataExpanded(*tempDataConst);      DataAbstract* temp=new DataExpanded(*tempDataConst);
900  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
901  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
902    } else if (isTagged()) {    } else if (isTagged()) {
903      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
904      DataAbstract* temp=new DataExpanded(*tempDataTag);      DataAbstract* temp=new DataExpanded(*tempDataTag);
905  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
906  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
907    } else if (isExpanded()) {    } else if (isExpanded()) {
908      //      //
909      // do nothing      // do nothing
# Line 656  Data::tag() Line 923  Data::tag()
923    if (isConstant()) {    if (isConstant()) {
924      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
925      DataAbstract* temp=new DataTagged(*tempDataConst);      DataAbstract* temp=new DataTagged(*tempDataConst);
926  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
927  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
928    } else if (isTagged()) {    } else if (isTagged()) {
929      // do nothing      // do nothing
930    } else if (isExpanded()) {    } else if (isExpanded()) {
# Line 671  Data::tag() Line 937  Data::tag()
937       {       {
938      throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");      throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");
939       }       }
940       m_data=res;      //      m_data=res;
941         set_m_data(res);
942       tag();       tag();
943    } else {    } else {
944      throw DataException("Error - Tagging not implemented for this Data type.");      throw DataException("Error - Tagging not implemented for this Data type.");
# Line 683  Data::resolve() Line 950  Data::resolve()
950  {  {
951    if (isLazy())    if (isLazy())
952    {    {
953       m_data=m_data->resolve();  //      m_data=m_data->resolve();
954        set_m_data(m_data->resolve());
955    }    }
956  }  }
957    
958    void
959    Data::requireWrite()
960    {
961      resolve();
962      exclusiveWrite();
963    }
964    
965  Data  Data
966  Data::oneOver() const  Data::oneOver() const
967  {  {
968    if (isLazy())    MAKELAZYOP(RECIP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),RECIP);  
     return Data(c);  
   }  
969    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
970  }  }
971    
972  Data  Data
973  Data::wherePositive() const  Data::wherePositive() const
974  {  {
975    if (isLazy())    MAKELAZYOP(GZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GZ);  
     return Data(c);  
   }  
976    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
977  }  }
978    
979  Data  Data
980  Data::whereNegative() const  Data::whereNegative() const
981  {  {
982    if (isLazy())    MAKELAZYOP(LZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LZ);  
     return Data(c);  
   }  
983    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
984  }  }
985    
986  Data  Data
987  Data::whereNonNegative() const  Data::whereNonNegative() const
988  {  {
989    if (isLazy())    MAKELAZYOP(GEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GEZ);  
     return Data(c);  
   }  
990    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
991  }  }
992    
993  Data  Data
994  Data::whereNonPositive() const  Data::whereNonPositive() const
995  {  {
996    if (isLazy())    MAKELAZYOP(LEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LEZ);  
     return Data(c);  
   }  
997    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
998  }  }
999    
1000  Data  Data
1001  Data::whereZero(double tol) const  Data::whereZero(double tol) const
1002  {  {
1003    Data dataAbs=abs();  //   Data dataAbs=abs();
1004    return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
1005       MAKELAZYOPOFF(EZ,tol)
1006       return C_TensorUnaryOperation(*this, bind2nd(AbsLTE(),tol));
1007    
1008  }  }
1009    
1010  Data  Data
1011  Data::whereNonZero(double tol) const  Data::whereNonZero(double tol) const
1012  {  {
1013    Data dataAbs=abs();  //   Data dataAbs=abs();
1014    return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
1015      MAKELAZYOPOFF(NEZ,tol)
1016      return C_TensorUnaryOperation(*this, bind2nd(AbsGT(),tol));
1017    
1018  }  }
1019    
1020  Data  Data
# Line 767  bool Line 1027  bool
1027  Data::probeInterpolation(const FunctionSpace& functionspace) const  Data::probeInterpolation(const FunctionSpace& functionspace) const
1028  {  {
1029    return getFunctionSpace().probeInterpolation(functionspace);    return getFunctionSpace().probeInterpolation(functionspace);
 //   if (getFunctionSpace()==functionspace) {  
 //     return true;  
 //   } else {  
 //     const_Domain_ptr domain=getDomain();  
 //     if  (*domain==*functionspace.getDomain()) {  
 //       return domain->probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());  
 //     } else {  
 //       return domain->probeInterpolationACross(getFunctionSpace().getTypeCode(),*(functionspace.getDomain()),functionspace.getTypeCode());  
 //     }  
 //   }  
1030  }  }
1031    
1032  Data  Data
# Line 813  Data::getDataPointSize() const Line 1063  Data::getDataPointSize() const
1063    return m_data->getNoValues();    return m_data->getNoValues();
1064  }  }
1065    
1066    
1067  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1068  Data::getLength() const  Data::getLength() const
1069  {  {
1070    return m_data->getLength();    return m_data->getLength();
1071  }  }
1072    
 const  
 boost::python::numeric::array  
 Data:: getValueOfDataPoint(int dataPointNo)  
 {  
   int i, j, k, l;  
   
   FORCERESOLVE;  
   
   //  
   // determine the rank and shape of each data point  
   int dataPointRank = getDataPointRank();  
   const DataTypes::ShapeType& dataPointShape = getDataPointShape();  
   
   //  
   // create the numeric array to be returned  
   boost::python::numeric::array numArray(0.0);  
1073    
1074    //  // There is no parallelism here ... elements need to be added in the correct order.
1075    // the shape of the returned numeric array will be the same  //   If we could presize the list and then fill in the elements it might work
1076    // as that of the data point  //   This would need setting elements to be threadsafe.
1077    int arrayRank = dataPointRank;  //   Having mulitple C threads calling into one interpreter is aparently a no-no.
1078    const DataTypes::ShapeType& arrayShape = dataPointShape;  const boost::python::object
1079    Data::toListOfTuples(bool scalarastuple)
1080    {
1081        using namespace boost::python;
1082        using boost::python::list;
1083        if (get_MPISize()>1)
1084        {
1085            throw DataException("::toListOfTuples is not available for MPI with more than one process.");
1086        }
1087        unsigned int rank=getDataPointRank();
1088        unsigned int size=getDataPointSize();
1089    
1090    //      int npoints=getNumDataPoints();
1091    // resize the numeric array to the shape just calculated      expand();           // This will also resolve if required
1092    if (arrayRank==0) {      const DataTypes::ValueType& vec=getReady()->getVectorRO();
1093      numArray.resize(1);      boost::python::list temp;
1094    }      temp.append(object());
1095    if (arrayRank==1) {      boost::python::list res(temp*npoints);// presize the list by the "[None] * npoints"  trick
1096      numArray.resize(arrayShape[0]);      if (rank==0)
1097    }      {
1098    if (arrayRank==2) {          long count;
1099      numArray.resize(arrayShape[0],arrayShape[1]);          if (scalarastuple)
1100    }          {
1101    if (arrayRank==3) {              for (count=0;count<npoints;++count)
1102      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);              {
1103    }          res[count]=make_tuple(vec[count]);
1104    if (arrayRank==4) {              }
1105      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);          }
1106    }          else
1107            {
1108                for (count=0;count<npoints;++count)
1109                {
1110                    res[count]=vec[count];
1111                }
1112            }
1113        }
1114        else if (rank==1)
1115        {
1116            size_t count;
1117            size_t offset=0;
1118            for (count=0;count<npoints;++count,offset+=size)
1119            {
1120                res[count]=pointToTuple1(getDataPointShape(), vec, offset);
1121            }
1122        }
1123        else if (rank==2)
1124        {
1125            size_t count;
1126            size_t offset=0;
1127            for (count=0;count<npoints;++count,offset+=size)
1128            {
1129            res[count]=pointToTuple2(getDataPointShape(), vec, offset);
1130            }
1131        }
1132        else if (rank==3)
1133        {
1134            size_t count;
1135            size_t offset=0;
1136            for (count=0;count<npoints;++count,offset+=size)
1137            {
1138                res[count]=pointToTuple3(getDataPointShape(), vec, offset);
1139            }
1140        }
1141        else if (rank==4)
1142        {
1143            size_t count;
1144            size_t offset=0;
1145            for (count=0;count<npoints;++count,offset+=size)
1146            {
1147                res[count]=pointToTuple4(getDataPointShape(), vec, offset);
1148            }
1149        }
1150        else
1151        {
1152            throw DataException("Unknown rank in ::toListOfTuples()");
1153        }
1154        return res;
1155    }
1156    
1157    if (getNumDataPointsPerSample()>0) {  const boost::python::object
1158    Data::getValueOfDataPointAsTuple(int dataPointNo)
1159    {
1160       forceResolve();
1161       if (getNumDataPointsPerSample()>0) {
1162         int sampleNo = dataPointNo/getNumDataPointsPerSample();         int sampleNo = dataPointNo/getNumDataPointsPerSample();
1163         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1164         //         //
1165         // Check a valid sample number has been supplied         // Check a valid sample number has been supplied
1166         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1167             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");             throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.");
1168         }         }
1169    
1170         //         //
1171         // Check a valid data point number has been supplied         // Check a valid data point number has been supplied
1172         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1173             throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");             throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid dataPointNoInSample.");
1174         }         }
1175         // TODO: global error handling         // TODO: global error handling
        // create a view of the data if it is stored locally  
 //       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);  
1176         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1177           return pointToTuple(getDataPointShape(),&(getDataAtOffsetRO(offset)));
1178      }
1179         switch( dataPointRank ){    else
1180              case 0 :    {
1181                  numArray[0] = getDataAtOffset(offset);      // The pre-numpy method would return an empty array of the given shape
1182                  break;      // I'm going to throw an exception because if we have zero points per sample we have problems
1183              case 1 :      throw DataException("Error - need at least 1 datapoint per sample.");
                 for( i=0; i<dataPointShape[0]; i++ )  
                     numArray[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));  
                 break;  
             case 2 :  
                 for( i=0; i<dataPointShape[0]; i++ )  
                     for( j=0; j<dataPointShape[1]; j++)  
                         numArray[make_tuple(i,j)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));  
                 break;  
             case 3 :  
                 for( i=0; i<dataPointShape[0]; i++ )  
                     for( j=0; j<dataPointShape[1]; j++ )  
                         for( k=0; k<dataPointShape[2]; k++)  
                             numArray[make_tuple(i,j,k)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));  
                 break;  
             case 4 :  
                 for( i=0; i<dataPointShape[0]; i++ )  
                     for( j=0; j<dataPointShape[1]; j++ )  
                         for( k=0; k<dataPointShape[2]; k++ )  
                             for( l=0; l<dataPointShape[3]; l++)  
                                 numArray[make_tuple(i,j,k,l)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));  
                 break;  
     }  
1184    }    }
   //  
   // return the array  
   return numArray;  
1185    
1186  }  }
1187    
1188    
1189  void  void
1190  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
1191  {  {
1192      // this will throw if the value cannot be represented      // this will throw if the value cannot be represented
1193      boost::python::numeric::array num_array(py_object);      setValueOfDataPointToArray(dataPointNo,py_object);
     setValueOfDataPointToArray(dataPointNo,num_array);  
1194  }  }
1195    
1196  void  void
1197  Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)  Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::object& obj)
1198  {  {
1199    if (isProtected()) {    if (isProtected()) {
1200          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1201    }    }
1202    FORCERESOLVE;    forceResolve();
1203    
1204      WrappedArray w(obj);
1205    //    //
1206    // check rank    // check rank
1207    if (static_cast<unsigned int>(num_array.getrank())<getDataPointRank())    if (static_cast<unsigned int>(w.getRank())<getDataPointRank())
1208        throw DataException("Rank of numarray does not match Data object rank");        throw DataException("Rank of array does not match Data object rank");
1209    
1210    //    //
1211    // check shape of num_array    // check shape of array
1212    for (unsigned int i=0; i<getDataPointRank(); i++) {    for (unsigned int i=0; i<getDataPointRank(); i++) {
1213      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])      if (w.getShape()[i]!=getDataPointShape()[i])
1214         throw DataException("Shape of numarray does not match Data object rank");         throw DataException("Shape of array does not match Data object rank");
1215    }    }
1216    //    //
1217    // make sure data is expanded:    // make sure data is expanded:
# Line 949  Data::setValueOfDataPointToArray(int dat Line 1222  Data::setValueOfDataPointToArray(int dat
1222    if (getNumDataPointsPerSample()>0) {    if (getNumDataPointsPerSample()>0) {
1223         int sampleNo = dataPointNo/getNumDataPointsPerSample();         int sampleNo = dataPointNo/getNumDataPointsPerSample();
1224         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1225         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,w);
1226    } else {    } else {
1227         m_data->copyToDataPoint(-1, 0,num_array);         m_data->copyToDataPoint(-1, 0,w);
1228    }    }
1229  }  }
1230    
# Line 963  Data::setValueOfDataPoint(int dataPointN Line 1236  Data::setValueOfDataPoint(int dataPointN
1236    }    }
1237    //    //
1238    // make sure data is expanded:    // make sure data is expanded:
1239    FORCERESOLVE;    forceResolve();
1240    if (!isExpanded()) {    if (!isExpanded()) {
1241      expand();      expand();
1242    }    }
# Line 977  Data::setValueOfDataPoint(int dataPointN Line 1250  Data::setValueOfDataPoint(int dataPointN
1250  }  }
1251    
1252  const  const
1253  boost::python::numeric::array  boost::python::object
1254  Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)  Data::getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo)
1255  {  {
1256    size_t length=0;    // This could be lazier than it is now
1257    int i, j, k, l, pos;    forceResolve();
   FORCERESOLVE;  
   //  
   // determine the rank and shape of each data point  
   int dataPointRank = getDataPointRank();  
   const DataTypes::ShapeType& dataPointShape = getDataPointShape();  
1258    
1259    //    // copy datapoint into a buffer
1260    // create the numeric array to be returned    // broadcast buffer to all nodes
1261    boost::python::numeric::array numArray(0.0);    // convert buffer to tuple
1262      // return tuple
1263    
1264    //    const DataTypes::ShapeType& dataPointShape = getDataPointShape();
1265    // the shape of the returned numeric array will be the same    size_t length=DataTypes::noValues(dataPointShape);
   // as that of the data point  
   int arrayRank = dataPointRank;  
   const DataTypes::ShapeType& arrayShape = dataPointShape;  
   
   //  
   // resize the numeric array to the shape just calculated  
   if (arrayRank==0) {  
     numArray.resize(1);  
   }  
   if (arrayRank==1) {  
     numArray.resize(arrayShape[0]);  
   }  
   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]);  
   }  
1266    
1267    // added for the MPI communication    // added for the MPI communication
   length=1;  
   for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];  
1268    double *tmpData = new double[length];    double *tmpData = new double[length];
1269    
1270    //    // updated for the MPI case
1271    // load the values for the data point into the numeric array.    if( get_MPIRank()==procNo ){
1272          if (getNumDataPointsPerSample()>0) {
1273      // updated for the MPI case      int sampleNo = dataPointNo/getNumDataPointsPerSample();
1274      if( get_MPIRank()==procNo ){      int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1275               if (getNumDataPointsPerSample()>0) {      //
1276                  int sampleNo = dataPointNo/getNumDataPointsPerSample();      // Check a valid sample number has been supplied
1277                  int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();      if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1278                  //          throw DataException("Error - Data::getValueOfGlobalDataPointAsTuple: invalid sampleNo.");
                 // Check a valid sample number has been supplied  
                 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {  
                   throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
                 }  
   
                 //  
                 // Check a valid data point number has been supplied  
                 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {  
                   throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");  
                 }  
                 // TODO: global error handling  
         // create a view of the data if it is stored locally  
         //DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);  
         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);  
   
         // pack the data from the view into tmpData for MPI communication  
         pos=0;  
         switch( dataPointRank ){  
             case 0 :  
                 tmpData[0] = getDataAtOffset(offset);  
                 break;  
             case 1 :  
                 for( i=0; i<dataPointShape[0]; i++ )  
                     tmpData[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));  
                 break;  
             case 2 :  
                 for( i=0; i<dataPointShape[0]; i++ )  
                     for( j=0; j<dataPointShape[1]; j++, pos++ )  
                         tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));  
                 break;  
             case 3 :  
                 for( i=0; i<dataPointShape[0]; i++ )  
                     for( j=0; j<dataPointShape[1]; j++ )  
                         for( k=0; k<dataPointShape[2]; k++, pos++ )  
                             tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));  
                 break;  
             case 4 :  
                 for( i=0; i<dataPointShape[0]; i++ )  
                     for( j=0; j<dataPointShape[1]; j++ )  
                         for( k=0; k<dataPointShape[2]; k++ )  
                             for( l=0; l<dataPointShape[3]; l++, pos++ )  
                                 tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));  
                 break;  
         }  
             }  
1279      }      }
1280          #ifdef PASO_MPI  
1281          // broadcast the data to all other processes      //
1282      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );      // Check a valid data point number has been supplied
1283          #endif      if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1284            throw DataException("Error - Data::getValueOfGlobalDataPointAsTuple: invalid dataPointNoInSample.");
     // unpack the data  
     switch( dataPointRank ){  
         case 0 :  
             numArray[0]=tmpData[0];  
             break;  
         case 1 :  
             for( i=0; i<dataPointShape[0]; i++ )  
                 numArray[i]=tmpData[i];  
             break;  
         case 2 :  
             for( i=0; i<dataPointShape[0]; i++ )  
                 for( j=0; j<dataPointShape[1]; j++ )  
                    numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];  
             break;  
         case 3 :  
             for( i=0; i<dataPointShape[0]; i++ )  
                 for( j=0; j<dataPointShape[1]; j++ )  
                     for( k=0; k<dataPointShape[2]; k++ )  
                         numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];  
             break;  
         case 4 :  
             for( i=0; i<dataPointShape[0]; i++ )  
                 for( j=0; j<dataPointShape[1]; j++ )  
                     for( k=0; k<dataPointShape[2]; k++ )  
                         for( l=0; l<dataPointShape[3]; l++ )  
                                 numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];  
             break;  
1285      }      }
1286        // TODO: global error handling
1287        DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1288    
1289        memcpy(tmpData,&(getDataAtOffsetRO(offset)),length*sizeof(double));
1290         }
1291      }
1292    #ifdef PASO_MPI
1293      // broadcast the data to all other processes
1294      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1295    #endif
1296    
1297      delete [] tmpData;    boost::python::tuple t=pointToTuple(dataPointShape,tmpData);
1298      delete [] tmpData;
1299    //    //
1300    // return the loaded array    // return the loaded array
1301    return numArray;    return t;
1302    
1303  }  }
1304    
1305    
1306  boost::python::numeric::array  boost::python::object
1307  Data::integrate_const() const  Data::integrateToTuple_const() const
1308  {  {
1309    if (isLazy())    if (isLazy())
1310    {    {
# Line 1127  Data::integrate_const() const Line 1313  Data::integrate_const() const
1313    return integrateWorker();    return integrateWorker();
1314  }  }
1315    
1316  boost::python::numeric::array  boost::python::object
1317  Data::integrate()  Data::integrateToTuple()
1318  {  {
1319    if (isLazy())    if (isLazy())
1320    {    {
1321      expand();      expand();
1322    }    }
1323    return integrateWorker();    return integrateWorker();
 }  
   
1324    
1325    }
1326    
1327  boost::python::numeric::array  boost::python::object
1328  Data::integrateWorker() const  Data::integrateWorker() const
1329  {  {
   int index;  
   int rank = getDataPointRank();  
1330    DataTypes::ShapeType shape = getDataPointShape();    DataTypes::ShapeType shape = getDataPointShape();
1331    int dataPointSize = getDataPointSize();    int dataPointSize = getDataPointSize();
1332    
# Line 1151  Data::integrateWorker() const Line 1334  Data::integrateWorker() const
1334    // calculate the integral values    // calculate the integral values
1335    vector<double> integrals(dataPointSize);    vector<double> integrals(dataPointSize);
1336    vector<double> integrals_local(dataPointSize);    vector<double> integrals_local(dataPointSize);
1337      const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1338      if (dom==0)
1339      {            
1340        throw DataException("Can not integrate over non-continuous domains.");
1341      }
1342  #ifdef PASO_MPI  #ifdef PASO_MPI
1343    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals_local,*this);    dom->setToIntegrals(integrals_local,*this);
1344    // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory    // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1345    double *tmp = new double[dataPointSize];    double *tmp = new double[dataPointSize];
1346    double *tmp_local = new double[dataPointSize];    double *tmp_local = new double[dataPointSize];
1347    for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }    for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1348    MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );    MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1349    for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }    for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1350      tuple result=pointToTuple(shape,tmp);
1351    delete[] tmp;    delete[] tmp;
1352    delete[] tmp_local;    delete[] tmp_local;
1353  #else  #else
1354    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);    dom->setToIntegrals(integrals,*this);
1355    /*  double *tmp = new double[dataPointSize];
1356      for (int i=0; i<dataPointSize; i++) { tmp[i]=integrals[i]; }*/
1357      tuple result=pointToTuple(shape,integrals);
1358    //   delete tmp;
1359  #endif  #endif
1360    
   //  
   // create the numeric array to be returned  
   // and load the array with the integral values  
   boost::python::numeric::array bp_array(1.0);  
   if (rank==0) {  
     bp_array.resize(1);  
     index = 0;  
     bp_array[0] = integrals[index];  
   }  
   if (rank==1) {  
     bp_array.resize(shape[0]);  
     for (int i=0; i<shape[0]; i++) {  
       index = i;  
       bp_array[i] = integrals[index];  
     }  
   }  
   if (rank==2) {  
        bp_array.resize(shape[0],shape[1]);  
        for (int i=0; i<shape[0]; i++) {  
          for (int j=0; j<shape[1]; j++) {  
            index = i + shape[0] * j;  
            bp_array[make_tuple(i,j)] = integrals[index];  
          }  
        }  
   }  
   if (rank==3) {  
     bp_array.resize(shape[0],shape[1],shape[2]);  
     for (int i=0; i<shape[0]; i++) {  
       for (int j=0; j<shape[1]; j++) {  
         for (int k=0; k<shape[2]; k++) {  
           index = i + shape[0] * ( j + shape[1] * k );  
           bp_array[make_tuple(i,j,k)] = integrals[index];  
         }  
       }  
     }  
   }  
   if (rank==4) {  
     bp_array.resize(shape[0],shape[1],shape[2],shape[3]);  
     for (int i=0; i<shape[0]; i++) {  
       for (int j=0; j<shape[1]; j++) {  
         for (int k=0; k<shape[2]; k++) {  
           for (int l=0; l<shape[3]; l++) {  
             index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );  
             bp_array[make_tuple(i,j,k,l)] = integrals[index];  
           }  
         }  
       }  
     }  
   }  
1361    
1362    //    return result;
   // return the loaded array  
   return bp_array;  
1363  }  }
1364    
1365  Data  Data
1366  Data::sin() const  Data::sin() const
1367  {  {
1368    if (isLazy())    MAKELAZYOP(SIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIN);  
     return Data(c);  
   }  
1369    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1370  }  }
1371    
1372  Data  Data
1373  Data::cos() const  Data::cos() const
1374  {  {
1375    if (isLazy())    MAKELAZYOP(COS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COS);  
     return Data(c);  
   }  
1376    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1377  }  }
1378    
1379  Data  Data
1380  Data::tan() const  Data::tan() const
1381  {  {
1382    if (isLazy())    MAKELAZYOP(TAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TAN);  
     return Data(c);  
   }  
1383    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1384  }  }
1385    
1386  Data  Data
1387  Data::asin() const  Data::asin() const
1388  {  {
1389    if (isLazy())    MAKELAZYOP(ASIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ASIN);  
     return Data(c);  
   }  
1390    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1391  }  }
1392    
1393  Data  Data
1394  Data::acos() const  Data::acos() const
1395  {  {
1396    if (isLazy())    MAKELAZYOP(ACOS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOS);  
     return Data(c);  
   }  
1397    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1398  }  }
1399    
# Line 1279  Data::acos() const Line 1401  Data::acos() const
1401  Data  Data
1402  Data::atan() const  Data::atan() const
1403  {  {
1404    if (isLazy())    MAKELAZYOP(ATAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ATAN);  
     return Data(c);  
   }  
1405    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1406  }  }
1407    
1408  Data  Data
1409  Data::sinh() const  Data::sinh() const
1410  {  {
1411    if (isLazy())    MAKELAZYOP(SINH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SINH);  
     return Data(c);  
   }  
1412    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1413  }  }
1414    
1415  Data  Data
1416  Data::cosh() const  Data::cosh() const
1417  {  {
1418    if (isLazy())    MAKELAZYOP(COSH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COSH);  
     return Data(c);  
   }  
1419    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1420  }  }
1421    
1422  Data  Data
1423  Data::tanh() const  Data::tanh() const
1424  {  {
1425    if (isLazy())    MAKELAZYOP(TANH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TANH);  
     return Data(c);  
   }  
1426    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1427  }  }
1428    
# Line 1324  Data::tanh() const Line 1430  Data::tanh() const
1430  Data  Data
1431  Data::erf() const  Data::erf() const
1432  {  {
1433  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1434    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1435  #else  #else
1436    if (isLazy())    MAKELAZYOP(ERF)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ERF);  
     return Data(c);  
   }  
1437    return C_TensorUnaryOperation(*this, ::erf);    return C_TensorUnaryOperation(*this, ::erf);
1438  #endif  #endif
1439  }  }
# Line 1339  Data::erf() const Line 1441  Data::erf() const
1441  Data  Data
1442  Data::asinh() const  Data::asinh() const
1443  {  {
1444    if (isLazy())    MAKELAZYOP(ASINH)
1445    {  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
     DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);  
     return Data(c);  
   }  
 #ifdef _WIN32  
1446    return C_TensorUnaryOperation(*this, escript::asinh_substitute);    return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1447  #else  #else
1448    return C_TensorUnaryOperation(*this, ::asinh);    return C_TensorUnaryOperation(*this, ::asinh);
# Line 1354  Data::asinh() const Line 1452  Data::asinh() const
1452  Data  Data
1453  Data::acosh() const  Data::acosh() const
1454  {  {
1455    if (isLazy())    MAKELAZYOP(ACOSH)
1456    {  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);  
     return Data(c);  
   }  
 #ifdef _WIN32  
1457    return C_TensorUnaryOperation(*this, escript::acosh_substitute);    return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1458  #else  #else
1459    return C_TensorUnaryOperation(*this, ::acosh);    return C_TensorUnaryOperation(*this, ::acosh);
# Line 1369  Data::acosh() const Line 1463  Data::acosh() const
1463  Data  Data
1464  Data::atanh() const  Data::atanh() const
1465  {  {
1466    if (isLazy())    MAKELAZYOP(ATANH)
1467    {  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
     DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);  
     return Data(c);  
   }  
 #ifdef _WIN32  
1468    return C_TensorUnaryOperation(*this, escript::atanh_substitute);    return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1469  #else  #else
1470    return C_TensorUnaryOperation(*this, ::atanh);    return C_TensorUnaryOperation(*this, ::atanh);
# Line 1383  Data::atanh() const Line 1473  Data::atanh() const
1473    
1474  Data  Data
1475  Data::log10() const  Data::log10() const
1476  {  if (isLazy())  {
1477    {    MAKELAZYOP(LOG10)
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG10);  
     return Data(c);  
   }  
1478    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1479  }  }
1480    
1481  Data  Data
1482  Data::log() const  Data::log() const
1483  {  {
1484    if (isLazy())    MAKELAZYOP(LOG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG);  
     return Data(c);  
   }  
1485    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1486  }  }
1487    
1488  Data  Data
1489  Data::sign() const  Data::sign() const
1490  {  {
1491    if (isLazy())    MAKELAZYOP(SIGN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIGN);  
     return Data(c);  
   }  
1492    return C_TensorUnaryOperation(*this, escript::fsign);    return C_TensorUnaryOperation(*this, escript::fsign);
1493  }  }
1494    
1495  Data  Data
1496  Data::abs() const  Data::abs() const
1497  {  {
1498    if (isLazy())    MAKELAZYOP(ABS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ABS);  
     return Data(c);  
   }  
1499    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1500  }  }
1501    
1502  Data  Data
1503  Data::neg() const  Data::neg() const
1504  {  {
1505    if (isLazy())    MAKELAZYOP(NEG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NEG);  
     return Data(c);  
   }  
1506    return C_TensorUnaryOperation(*this, negate<double>());    return C_TensorUnaryOperation(*this, negate<double>());
1507  }  }
1508    
# Line 1448  Data::pos() const Line 1519  Data::pos() const
1519    
1520  Data  Data
1521  Data::exp() const  Data::exp() const
1522  {    {
1523    if (isLazy())    MAKELAZYOP(EXP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),EXP);  
     return Data(c);  
   }  
1524    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1525  }  }
1526    
1527  Data  Data
1528  Data::sqrt() const  Data::sqrt() const
1529  {  {
1530    if (isLazy())    MAKELAZYOP(SQRT)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SQRT);  
     return Data(c);  
   }  
1531    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1532  }  }
1533    
# Line 1483  Data::Lsup() Line 1546  Data::Lsup()
1546  {  {
1547     if (isLazy())     if (isLazy())
1548     {     {
1549      expand();      resolve();
1550     }     }
1551     return LsupWorker();     return LsupWorker();
1552  }  }
# Line 1503  Data::sup() Line 1566  Data::sup()
1566  {  {
1567     if (isLazy())     if (isLazy())
1568     {     {
1569      expand();      resolve();
1570     }     }
1571     return supWorker();     return supWorker();
1572  }  }
# Line 1523  Data::inf() Line 1586  Data::inf()
1586  {  {
1587     if (isLazy())     if (isLazy())
1588     {     {
1589      expand();      resolve();
1590     }     }
1591     return infWorker();     return infWorker();
1592  }  }
# Line 1633  Data::swapaxes(const int axis0, const in Line 1696  Data::swapaxes(const int axis0, const in
1696       if (axis0 == axis1) {       if (axis0 == axis1) {
1697           throw DataException("Error - Data::swapaxes: axis indices must be different.");           throw DataException("Error - Data::swapaxes: axis indices must be different.");
1698       }       }
1699       if (axis0 > axis1) {       MAKELAZYOP2(SWAP,axis0,axis1)
1700           axis0_tmp=axis1;       if (axis0 > axis1)
1701           axis1_tmp=axis0;       {
1702       } else {      axis0_tmp=axis1;
1703           axis0_tmp=axis0;      axis1_tmp=axis0;
1704           axis1_tmp=axis1;       }
1705         else
1706         {
1707        axis0_tmp=axis0;
1708        axis1_tmp=axis1;
1709       }       }
1710       for (int i=0; i<rank; i++) {       for (int i=0; i<rank; i++)
1711         if (i == axis0_tmp) {       {
1712            ev_shape.push_back(s[axis1_tmp]);      if (i == axis0_tmp)
1713         } else if (i == axis1_tmp) {      {
1714            ev_shape.push_back(s[axis0_tmp]);          ev_shape.push_back(s[axis1_tmp]);
1715         } else {      }
1716            ev_shape.push_back(s[i]);      else if (i == axis1_tmp)
1717         }      {
1718            ev_shape.push_back(s[axis0_tmp]);
1719        }
1720        else
1721        {
1722            ev_shape.push_back(s[i]);
1723        }
1724       }       }
1725       Data ev(0.,ev_shape,getFunctionSpace());       Data ev(0.,ev_shape,getFunctionSpace());
1726       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1727       m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);       m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1728       return ev;       return ev;
   
1729  }  }
1730    
1731  Data  Data
# Line 1672  Data::symmetric() const Line 1744  Data::symmetric() const
1744       else {       else {
1745          throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");          throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1746       }       }
1747       if (isLazy())       MAKELAZYOP(SYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SYM);  
     return Data(c);  
      }  
1748       Data ev(0.,getDataPointShape(),getFunctionSpace());       Data ev(0.,getDataPointShape(),getFunctionSpace());
1749       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1750       m_data->symmetric(ev.m_data.get());       m_data->symmetric(ev.m_data.get());
# Line 1686  Data::symmetric() const Line 1754  Data::symmetric() const
1754  Data  Data
1755  Data::nonsymmetric() const  Data::nonsymmetric() const
1756  {  {
1757       if (isLazy())       MAKELAZYOP(NSYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NSYM);  
     return Data(c);  
      }  
1758       // check input       // check input
1759       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1760       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1722  Data::nonsymmetric() const Line 1786  Data::nonsymmetric() const
1786       }       }
1787  }  }
1788    
   
 // Doing a lazy version of this would require some thought.  
 // First it needs a parameter (which DataLazy doesn't support at the moment).  
 // (secondly although it does not apply to trace) we can't handle operations which return  
 // multiple results (like eigenvectors_values) or return values of different shapes to their input  
 // (like eigenvalues).  
1789  Data  Data
1790  Data::trace(int axis_offset) const  Data::trace(int axis_offset) const
1791  {  {    
1792       if (isLazy())       MAKELAZYOPOFF(TRACE,axis_offset)
1793         if ((axis_offset<0) || (axis_offset>getDataPointRank()))
1794       {       {
1795      Data temp(*this);   // to get around the fact that you can't resolve a const Data      throw DataException("Error - Data::trace, axis_offset must be between 0 and rank-2 inclusive.");
     temp.resolve();  
     return temp.trace(axis_offset);  
1796       }       }
1797       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1798       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1787  Data::trace(int axis_offset) const Line 1844  Data::trace(int axis_offset) const
1844  Data  Data
1845  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1846  {      {    
1847       if (isLazy())       MAKELAZYOPOFF(TRANS,axis_offset)
      {  
     Data temp(*this);   // to get around the fact that you can't resolve a const Data  
     temp.resolve();  
     return temp.transpose(axis_offset);  
      }  
1848       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1849       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
1850       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
# Line 1894  Data::calc_minGlobalDataPoint(int& ProcN Line 1946  Data::calc_minGlobalDataPoint(int& ProcN
1946    double next,local_min;    double next,local_min;
1947    int local_lowi=0,local_lowj=0;        int local_lowi=0,local_lowj=0;    
1948    
1949    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)    #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
1950    {    {
1951      local_min=min;      local_min=min;
1952      #pragma omp for private(i,j) schedule(static)      #pragma omp for private(i,j) schedule(static)
1953      for (i=0; i<numSamples; i++) {      for (i=0; i<numSamples; i++) {
1954        for (j=0; j<numDPPSample; j++) {        for (j=0; j<numDPPSample; j++) {
1955          next=temp.getDataAtOffset(temp.getDataOffset(i,j));          next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
1956          if (next<local_min) {          if (next<local_min) {
1957            local_min=next;            local_min=next;
1958            local_lowi=i;            local_lowi=i;
# Line 1909  Data::calc_minGlobalDataPoint(int& ProcN Line 1961  Data::calc_minGlobalDataPoint(int& ProcN
1961        }        }
1962      }      }
1963      #pragma omp critical      #pragma omp critical
1964      if (local_min<min) {      if (local_min<min) {    // If we found a smaller value than our sentinel
1965        min=local_min;        min=local_min;
1966        lowi=local_lowi;        lowi=local_lowi;
1967        lowj=local_lowj;        lowj=local_lowj;
# Line 1917  Data::calc_minGlobalDataPoint(int& ProcN Line 1969  Data::calc_minGlobalDataPoint(int& ProcN
1969    }    }
1970    
1971  #ifdef PASO_MPI  #ifdef PASO_MPI
1972      // determine the processor on which the minimum occurs    // determine the processor on which the minimum occurs
1973      next = temp.getDataPoint(lowi,lowj);    next = temp.getDataPointRO(lowi,lowj);
1974      int lowProc = 0;    int lowProc = 0;
1975      double *globalMins = new double[get_MPISize()+1];    double *globalMins = new double[get_MPISize()+1];
1976      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );    int error;
1977      error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1978      if( get_MPIRank()==0 ){  
1979          next = globalMins[lowProc];    if( get_MPIRank()==0 ){
1980          for( i=1; i<get_MPISize(); i++ )      next = globalMins[lowProc];
1981              if( next>globalMins[i] ){      for( i=1; i<get_MPISize(); i++ )
1982                  lowProc = i;          if( next>globalMins[i] ){
1983                  next = globalMins[i];              lowProc = i;
1984              }              next = globalMins[i];
1985            }
1986      }
1987      MPI_Bcast( &lowProc, 1, MPI_INT, 0, get_MPIComm() );
1988      DataPointNo = lowj + lowi * numDPPSample;
1989      MPI_Bcast(&DataPointNo, 1, MPI_INT, lowProc, get_MPIComm() );
1990      delete [] globalMins;
1991      ProcNo = lowProc;
1992    #else
1993      ProcNo = 0;
1994      DataPointNo = lowj + lowi * numDPPSample;
1995    #endif
1996    }
1997    
1998    
1999    const boost::python::tuple
2000    Data::maxGlobalDataPoint() const
2001    {
2002      int DataPointNo;
2003      int ProcNo;
2004      calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2005      return make_tuple(ProcNo,DataPointNo);
2006    }
2007    
2008    void
2009    Data::calc_maxGlobalDataPoint(int& ProcNo,
2010                            int& DataPointNo) const
2011    {
2012      if (isLazy())
2013      {
2014        Data temp(*this);   // to get around the fact that you can't resolve a const Data
2015        temp.resolve();
2016        return temp.calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2017      }
2018      int i,j;
2019      int highi=0,highj=0;
2020    //-------------
2021      double max= -numeric_limits<double>::max();
2022    
2023      Data temp=maxval();
2024    
2025      int numSamples=temp.getNumSamples();
2026      int numDPPSample=temp.getNumDataPointsPerSample();
2027    
2028      double next,local_max;
2029      int local_highi=0,local_highj=0;  
2030    
2031      #pragma omp parallel firstprivate(local_highi,local_highj) private(next,local_max)
2032      {
2033        local_max=max;
2034        #pragma omp for private(i,j) schedule(static)
2035        for (i=0; i<numSamples; i++) {
2036          for (j=0; j<numDPPSample; j++) {
2037            next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2038            if (next>local_max) {
2039              local_max=next;
2040              local_highi=i;
2041              local_highj=j;
2042            }
2043          }
2044        }
2045        #pragma omp critical
2046        if (local_max>max) {    // If we found a larger value than our sentinel
2047          max=local_max;
2048          highi=local_highi;
2049          highj=local_highj;
2050        }
2051      }
2052    #ifdef PASO_MPI
2053      // determine the processor on which the maximum occurs
2054      next = temp.getDataPointRO(highi,highj);
2055      int highProc = 0;
2056      double *globalMaxs = new double[get_MPISize()+1];
2057      int error;
2058      error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMaxs, 1, MPI_DOUBLE, 0, get_MPIComm() );
2059      if( get_MPIRank()==0 ){
2060        next = globalMaxs[highProc];
2061        for( i=1; i<get_MPISize(); i++ )
2062        {
2063        if( next<globalMaxs[i] )
2064        {
2065            highProc = i;
2066            next = globalMaxs[i];
2067      }      }
2068      MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );      }
2069      }
2070      MPI_Bcast( &highProc, 1, MPI_INT, 0, get_MPIComm() );
2071      DataPointNo = highj + highi * numDPPSample;  
2072      MPI_Bcast(&DataPointNo, 1, MPI_INT, highProc, get_MPIComm() );
2073    
2074      delete [] globalMins;    delete [] globalMaxs;
2075      ProcNo = lowProc;    ProcNo = highProc;
2076  #else  #else
2077      ProcNo = 0;    ProcNo = 0;
2078      DataPointNo = highj + highi * numDPPSample;
2079  #endif  #endif
   DataPointNo = lowj + lowi * numDPPSample;  
2080  }  }
2081    
2082  void  void
# Line 1977  Data::saveVTK(std::string fileName) cons Line 2115  Data::saveVTK(std::string fileName) cons
2115    }    }
2116    boost::python::dict args;    boost::python::dict args;
2117    args["data"]=boost::python::object(this);    args["data"]=boost::python::object(this);
2118    getDomain()->saveVTK(fileName,args);    getDomain()->saveVTK(fileName,args,"","");
2119    return;    return;
2120  }  }
2121    
2122    
2123    
2124  Data&  Data&
2125  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
2126  {  {
2127    if (isProtected()) {    if (isProtected()) {
2128          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2129    }    }
2130    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,ADD)    // for lazy + is equivalent to +=
2131    {    exclusiveWrite();         // Since Lazy data does not modify its leaves we only need to worry here
2132      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),ADD); // for lazy + is equivalent to +=    binaryOp(right,plus<double>());
2133          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,plus<double>());  
     return (*this);  
   }  
2134  }  }
2135    
2136  Data&  Data&
# Line 2007  Data::operator+=(const boost::python::ob Line 2140  Data::operator+=(const boost::python::ob
2140          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2141    }    }
2142    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2143    if (isLazy())    (*this)+=tmp;
2144    {    return *this;
     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),ADD);   // for lazy + is equivalent to +=  
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(tmp,plus<double>());  
     return (*this);  
   }  
2145  }  }
2146    
2147  // Hmmm, operator= makes a deep copy but the copy constructor does not?  // Hmmm, operator= makes a deep copy but the copy constructor does not?
2148  Data&  Data&
2149  Data::operator=(const Data& other)  Data::operator=(const Data& other)
2150  {  {
2151    copy(other);    m_protected=false;        // since any changes should be caught by exclusiveWrite();
2152    //   m_data=other.m_data;
2153      set_m_data(other.m_data);
2154    return (*this);    return (*this);
2155  }  }
2156    
# Line 2034  Data::operator-=(const Data& right) Line 2160  Data::operator-=(const Data& right)
2160    if (isProtected()) {    if (isProtected()) {
2161          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2162    }    }
2163    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,SUB)
2164    {    exclusiveWrite();
2165      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),SUB); // for lazy - is equivalent to -=    binaryOp(right,minus<double>());
2166          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,minus<double>());  
     return (*this);  
   }  
2167  }  }
2168    
2169  Data&  Data&
# Line 2054  Data::operator-=(const boost::python::ob Line 2173  Data::operator-=(const boost::python::ob
2173          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2174    }    }
2175    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2176    if (isLazy())    (*this)-=tmp;
2177    {    return (*this);
     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),SUB);   // for lazy - is equivalent to -=  
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(tmp,minus<double>());  
     return (*this);  
   }  
2178  }  }
2179    
2180  Data&  Data&
# Line 2073  Data::operator*=(const Data& right) Line 2183  Data::operator*=(const Data& right)
2183    if (isProtected()) {    if (isProtected()) {
2184          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2185    }    }
2186    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,MUL)
2187    {    exclusiveWrite();
2188      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),MUL); // for lazy * is equivalent to *=    binaryOp(right,multiplies<double>());
2189          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,multiplies<double>());  
     return (*this);  
   }  
2190  }  }
2191    
2192  Data&  Data&
# Line 2093  Data::operator*=(const boost::python::ob Line 2196  Data::operator*=(const boost::python::ob
2196          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2197    }    }
2198    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2199    if (isLazy())    (*this)*=tmp;
2200    {    return (*this);
     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),MUL);   // for lazy * is equivalent to *=  
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(tmp,multiplies<double>());  
     return (*this);  
   }  
2201  }  }
2202    
2203  Data&  Data&
# Line 2112  Data::operator/=(const Data& right) Line 2206  Data::operator/=(const Data& right)
2206    if (isProtected()) {    if (isProtected()) {
2207          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2208    }    }
2209    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,DIV)
2210    {    exclusiveWrite();
2211      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),DIV); // for lazy / is equivalent to /=    binaryOp(right,divides<double>());
2212          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,divides<double>());  
     return (*this);  
   }  
2213  }  }
2214    
2215  Data&  Data&
# Line 2132  Data::operator/=(const boost::python::ob Line 2219  Data::operator/=(const boost::python::ob
2219          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2220    }    }
2221    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2222    if (isLazy())    (*this)/=tmp;
2223    {    return (*this);
     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),DIV);   // for lazy / is equivalent to /=  
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(tmp,divides<double>());  
     return (*this);  
   }  
2224  }  }
2225    
2226  Data  Data
# Line 2162  Data::powO(const boost::python::object& Line 2240  Data::powO(const boost::python::object&
2240  Data  Data
2241  Data::powD(const Data& right) const  Data::powD(const Data& right) const
2242  {  {
2243    if (isLazy() || right.isLazy())    MAKELAZYBIN(right,POW)
   {  
     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),POW);  
     return Data(c);  
   }  
2244    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2245  }  }
2246    
# Line 2175  Data::powD(const Data& right) const Line 2249  Data::powD(const Data& right) const
2249  Data  Data
2250  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
2251  {  {
2252    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,ADD)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);  
     return Data(c);  
   }  
2253    return C_TensorBinaryOperation(left, right, plus<double>());    return C_TensorBinaryOperation(left, right, plus<double>());
2254  }  }
2255    
# Line 2188  escript::operator+(const Data& left, con Line 2258  escript::operator+(const Data& left, con
2258  Data  Data
2259  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
2260  {  {
2261    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,SUB)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);  
     return Data(c);  
   }  
2262    return C_TensorBinaryOperation(left, right, minus<double>());    return C_TensorBinaryOperation(left, right, minus<double>());
2263  }  }
2264    
# Line 2201  escript::operator-(const Data& left, con Line 2267  escript::operator-(const Data& left, con
2267  Data  Data
2268  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
2269  {  {
2270    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,MUL)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);  
     return Data(c);  
   }  
2271    return C_TensorBinaryOperation(left, right, multiplies<double>());    return C_TensorBinaryOperation(left, right, multiplies<double>());
2272  }  }
2273    
# Line 2214  escript::operator*(const Data& left, con Line 2276  escript::operator*(const Data& left, con
2276  Data  Data
2277  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
2278  {  {
2279    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,DIV)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),DIV);  
     return Data(c);  
   }  
2280    return C_TensorBinaryOperation(left, right, divides<double>());    return C_TensorBinaryOperation(left, right, divides<double>());
2281  }  }
2282    
# Line 2227  escript::operator/(const Data& left, con Line 2285  escript::operator/(const Data& left, con
2285  Data  Data
2286  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
2287  {  {
2288    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2289    {    MAKELAZYBIN2(left,tmp,ADD)
2290      DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),ADD);    return left+tmp;
     return Data(c);  
   }  
   return left+Data(right,left.getFunctionSpace(),false);  
2291  }  }
2292    
2293  //  //
# Line 2240  escript::operator+(const Data& left, con Line 2295  escript::operator+(const Data& left, con
2295  Data  Data
2296  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
2297  {  {
2298    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2299    {    MAKELAZYBIN2(left,tmp,SUB)
2300      DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),SUB);    return left-tmp;
     return Data(c);  
   }  
   return left-Data(right,left.getFunctionSpace(),false);  
2301  }  }
2302    
2303  //  //
# Line 2253  escript::operator-(const Data& left, con Line 2305  escript::operator-(const Data& left, con
2305  Data  Data
2306  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
2307  {  {
2308    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2309    {    MAKELAZYBIN2(left,tmp,MUL)
2310      DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),MUL);    return left*tmp;
     return Data(c);  
   }  
   return left*Data(right,left.getFunctionSpace(),false);  
2311  }  }
2312    
2313  //  //
# Line 2266  escript::operator*(const Data& left, con Line 2315  escript::operator*(const Data& left, con
2315  Data  Data
2316  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
2317  {  {
2318    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2319    {    MAKELAZYBIN2(left,tmp,DIV)
2320      DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),DIV);    return left/tmp;
     return Data(c);  
   }  
   return left/Data(right,left.getFunctionSpace(),false);  
2321  }  }
2322    
2323  //  //
# Line 2279  escript::operator/(const Data& left, con Line 2325  escript::operator/(const Data& left, con
2325  Data  Data
2326  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
2327  {  {
2328    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2329    {    MAKELAZYBIN2(tmp,right,ADD)
2330      DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),ADD);    return tmp+right;
     return Data(c);  
   }  
   return Data(left,right.getFunctionSpace(),false)+right;  
2331  }  }
2332    
2333  //  //
# Line 2292  escript::operator+(const boost::python:: Line 2335  escript::operator+(const boost::python::
2335  Data  Data
2336  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
2337  {  {
2338    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2339    {    MAKELAZYBIN2(tmp,right,SUB)
2340      DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),SUB);    return tmp-right;
     return Data(c);  
   }  
   return Data(left,right.getFunctionSpace(),false)-right;  
2341  }  }
2342    
2343  //  //
# Line 2305  escript::operator-(const boost::python:: Line 2345  escript::operator-(const boost::python::
2345  Data  Data
2346  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
2347  {  {
2348    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2349    {    MAKELAZYBIN2(tmp,right,MUL)
2350      DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),MUL);    return tmp*right;
     return Data(c);  
   }  
   return Data(left,right.getFunctionSpace(),false)*right;  
2351  }  }
2352    
2353  //  //
# Line 2318  escript::operator*(const boost::python:: Line 2355  escript::operator*(const boost::python::
2355  Data  Data
2356  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
2357  {  {
2358    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2359    {    MAKELAZYBIN2(tmp,right,DIV)
2360      DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),DIV);    return tmp/right;
     return Data(c);  
   }  
   return Data(left,right.getFunctionSpace(),false)/right;  
2361  }  }
2362    
2363    
# Line 2364  void Line 2398  void
2398  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
2399                 const Data& value)                 const Data& value)
2400  {  {
 //  const DataArrayView& view=getPointDataView();  
   
2401    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2402    if (slice_region.size()!=getDataPointRank()) {    if (slice_region.size()!=getDataPointRank()) {
2403      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
2404    }    }
2405      exclusiveWrite();
2406    if (getFunctionSpace()!=value.getFunctionSpace()) {    if (getFunctionSpace()!=value.getFunctionSpace()) {
2407       setSlice(Data(value,getFunctionSpace()),slice_region);       setSlice(Data(value,getFunctionSpace()),slice_region);
2408    } else {    } else {
# Line 2384  Data::setSlice(const Data& value, Line 2417  Data::setSlice(const Data& value,
2417    if (isProtected()) {    if (isProtected()) {
2418          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2419    }    }
2420    FORCERESOLVE;    forceResolve();
2421  /*  if (isLazy())    exclusiveWrite();     // In case someone finds a way to call this without going through setItemD
   {  
     throw DataException("Error - setSlice not permitted on lazy data.");  
   }*/  
2422    Data tempValue(value);    Data tempValue(value);
2423    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
2424    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 2431  Data::typeMatchRight(const Data& right) Line 2461  Data::typeMatchRight(const Data& right)
2461    }    }
2462  }  }
2463    
2464    // The normal TaggedValue adds the tag if it is not already present
2465    // This form does not. It throws instead.
2466    // This is because the names are maintained by the domain and cannot be added
2467    // without knowing the tag number to map it to.
2468  void  void
2469  Data::setTaggedValueByName(std::string name,  Data::setTaggedValueByName(std::string name,
2470                             const boost::python::object& value)                             const boost::python::object& value)
2471  {  {
2472       if (getFunctionSpace().getDomain()->isValidTagName(name)) {       if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2473      FORCERESOLVE;      forceResolve();
2474        exclusiveWrite();
2475          int tagKey=getFunctionSpace().getDomain()->getTag(name);          int tagKey=getFunctionSpace().getDomain()->getTag(name);
2476          setTaggedValue(tagKey,value);          setTaggedValue(tagKey,value);
2477       }       }
2478         else
2479         {                  // The
2480        throw DataException("Error - unknown tag in setTaggedValueByName.");
2481         }
2482  }  }
2483    
2484  void  void
2485  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
2486                       const boost::python::object& value)                       const boost::python::object& value)
# Line 2450  Data::setTaggedValue(int tagKey, Line 2490  Data::setTaggedValue(int tagKey,
2490    }    }
2491    //    //
2492    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2493    FORCERESOLVE;    forceResolve();
2494      exclusiveWrite();
2495    if (isConstant()) tag();    if (isConstant()) tag();
2496    numeric::array asNumArray(value);    WrappedArray w(value);
   
   // extract the shape of the numarray  
   DataTypes::ShapeType tempShape;  
   for (int i=0; i < asNumArray.getrank(); i++) {  
     tempShape.push_back(extract<int>(asNumArray.getshape()[i]));  
   }  
2497    
2498    DataVector temp_data2;    DataVector temp_data2;
2499    temp_data2.copyFromNumArray(asNumArray);    temp_data2.copyFromArray(w,1);
2500    
2501    m_data->setTaggedValue(tagKey,tempShape, temp_data2);    m_data->setTaggedValue(tagKey,w.getShape(), temp_data2);
2502  }  }
2503    
2504    
# Line 2478  Data::setTaggedValueFromCPP(int tagKey, Line 2513  Data::setTaggedValueFromCPP(int tagKey,
2513    }    }
2514    //    //
2515    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2516    FORCERESOLVE;    forceResolve();
2517    if (isConstant()) tag();    if (isConstant()) tag();
2518      exclusiveWrite();
2519    //    //
2520    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
2521    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
# Line 2512  escript::C_GeneralTensorProduct(Data& ar Line 2548  escript::C_GeneralTensorProduct(Data& ar
2548    // SM is the product of the last axis_offset entries in arg_0.getShape().    // SM is the product of the last axis_offset entries in arg_0.getShape().
2549    
2550    // deal with any lazy data    // deal with any lazy data
2551    if (arg_0.isLazy()) {arg_0.resolve();}  //   if (arg_0.isLazy()) {arg_0.resolve();}
2552    if (arg_1.isLazy()) {arg_1.resolve();}  //   if (arg_1.isLazy()) {arg_1.resolve();}
2553      if (arg_0.isLazy() || arg_1.isLazy() || (AUTOLAZYON && (arg_0.isExpanded() || arg_1.isExpanded())))
2554      {
2555        DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2556        return Data(c);
2557      }
2558    
2559    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2560    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
# Line 2590  escript::C_GeneralTensorProduct(Data& ar Line 2631  escript::C_GeneralTensorProduct(Data& ar
2631       for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z       for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2632    }    }
2633    
2634      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2635      {
2636         ostringstream os;
2637         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2638         throw DataException(os.str());
2639      }
2640    
2641    // Declare output Data object    // Declare output Data object
2642    Data res;    Data res;
2643    
2644    if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {    if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2645      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2646      double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2647      double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2648      double *ptr_2 = &(res.getDataAtOffset(0));      double *ptr_2 = &(res.getDataAtOffsetRW(0));
2649      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2650    }    }
2651    else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {    else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
# Line 2618  escript::C_GeneralTensorProduct(Data& ar Line 2666  escript::C_GeneralTensorProduct(Data& ar
2666    
2667      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2668      int offset_0 = tmp_0->getPointOffset(0,0);      int offset_0 = tmp_0->getPointOffset(0,0);
2669      double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
     // Get the views  
 //     DataArrayView view_1 = tmp_1->getDefaultValue();  
 //     DataArrayView view_2 = tmp_2->getDefaultValue();  
 //     // Get the pointers to the actual data  
 //     double *ptr_1 = &((view_1.getData())[0]);  
 //     double *ptr_2 = &((view_2.getData())[0]);  
   
     double *ptr_1 = &(tmp_1->getDefaultValue(0));  
     double *ptr_2 = &(tmp_2->getDefaultValue(0));  
2670    
2671        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2672        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2673    
2674      // Compute an MVP for the default      // Compute an MVP for the default
2675      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
# Line 2637  escript::C_GeneralTensorProduct(Data& ar Line 2678  escript::C_GeneralTensorProduct(Data& ar
2678      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2679      for (i=lookup_1.begin();i!=lookup_1.end();i++) {      for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2680        tmp_2->addTag(i->first);        tmp_2->addTag(i->first);
 //       DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);  
 //       DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);  
 //       double *ptr_1 = &view_1.getData(0);  
 //       double *ptr_2 = &view_2.getData(0);  
2681    
2682        double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2683        double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2684            
2685        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2686      }      }
# Line 2667  escript::C_GeneralTensorProduct(Data& ar Line 2704  escript::C_GeneralTensorProduct(Data& ar
2704        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2705          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2706          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2707          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2708          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2709          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2710          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2711        }        }
2712      }      }
# Line 2693  escript::C_GeneralTensorProduct(Data& ar Line 2730  escript::C_GeneralTensorProduct(Data& ar
2730    
2731      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2732      int offset_1 = tmp_1->getPointOffset(0,0);      int offset_1 = tmp_1->getPointOffset(0,0);
2733      double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2734      // Get the views      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2735  //     DataArrayView view_0 = tmp_0->getDefaultValue();      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
 //     DataArrayView view_2 = tmp_2->getDefaultValue();  
 //     // Get the pointers to the actual data  
 //     double *ptr_0 = &((view_0.getData())[0]);  
 //     double *ptr_2 = &((view_2.getData())[0]);  
   
     double *ptr_0 = &(tmp_0->getDefaultValue(0));  
     double *ptr_2 = &(tmp_2->getDefaultValue(0));  
2736    
2737      // Compute an MVP for the default      // Compute an MVP for the default
2738      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
# Line 2710  escript::C_GeneralTensorProduct(Data& ar Line 2740  escript::C_GeneralTensorProduct(Data& ar
2740      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2741      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2742      for (i=lookup_0.begin();i!=lookup_0.end();i++) {      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
 //      tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());  
 //       DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);  
 //       DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);  
 //       double *ptr_0 = &view_0.getData(0);  
 //       double *ptr_2 = &view_2.getData(0);  
2743    
2744        tmp_2->addTag(i->first);        tmp_2->addTag(i->first);
2745        double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2746        double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2747        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2748      }      }
2749    
# Line 2739  escript::C_GeneralTensorProduct(Data& ar Line 2764  escript::C_GeneralTensorProduct(Data& ar
2764      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2765      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2766    
2767  //     // Get the views      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2768  //     DataArrayView view_0 = tmp_0->getDefaultValue();      const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2769  //     DataArrayView view_1 = tmp_1->getDefaultValue();      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
 //     DataArrayView view_2 = tmp_2->getDefaultValue();  
 //     // Get the pointers to the actual data  
 //     double *ptr_0 = &((view_0.getData())[0]);  
 //     double *ptr_1 = &((view_1.getData())[0]);  
 //     double *ptr_2 = &((view_2.getData())[0]);  
   
     double *ptr_0 = &(tmp_0->getDefaultValue(0));  
     double *ptr_1 = &(tmp_1->getDefaultValue(0));  
     double *ptr_2 = &(tmp_2->getDefaultValue(0));  
   
2770    
2771      // Compute an MVP for the default      // Compute an MVP for the default
2772      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
# Line 2768  escript::C_GeneralTensorProduct(Data& ar Line 2783  escript::C_GeneralTensorProduct(Data& ar
2783      // Compute an MVP for each tag      // Compute an MVP for each tag
2784      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2785      for (i=lookup_2.begin();i!=lookup_2.end();i++) {      for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2786  //       DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2787  //       DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2788  //       DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
 //       double *ptr_0 = &view_0.getData(0);  
 //       double *ptr_1 = &view_1.getData(0);  
 //       double *ptr_2 = &view_2.getData(0);  
   
       double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));  
       double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));  
       double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));  
2789    
2790        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2791      }      }
# Line 2799  escript::C_GeneralTensorProduct(Data& ar Line 2807  escript::C_GeneralTensorProduct(Data& ar
2807      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2808      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2809        int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0        int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2810        double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2811        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2812          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2813          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2814          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2815          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2816          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2817        }        }
2818      }      }
# Line 2828  escript::C_GeneralTensorProduct(Data& ar Line 2836  escript::C_GeneralTensorProduct(Data& ar
2836        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2837          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2838          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2839          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2840          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2841          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2842          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2843        }        }
2844      }      }
# Line 2853  escript::C_GeneralTensorProduct(Data& ar Line 2861  escript::C_GeneralTensorProduct(Data& ar
2861      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2862      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2863        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2864        double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2865        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2866          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2867          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2868          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2869          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2870          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2871        }        }
2872      }      }
# Line 2883  escript::C_GeneralTensorProduct(Data& ar Line 2891  escript::C_GeneralTensorProduct(Data& ar
2891          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2892          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2893          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2894          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2895          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2896          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2897          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2898        }        }
2899      }      }
# Line 2923  Data::borrowReadyPtr() const Line 2931  Data::borrowReadyPtr() const
2931  std::string  std::string
2932  Data::toString() const  Data::toString() const
2933  {  {
2934      if (!m_data->isEmpty() &&      if (!m_data->isEmpty() &&
2935      getNumDataPoints()*getDataPointSize()>escriptParams.getInt("TOO_MANY_LINES"))      !m_data->isLazy() &&
2936        getLength()>escriptParams.getInt("TOO_MANY_LINES"))
2937      {      {
2938      stringstream temp;      stringstream temp;
2939      temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();      temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();
# Line 2934  Data::toString() const Line 2943  Data::toString() const
2943  }  }
2944    
2945    
2946    // This method is not thread-safe
2947    DataTypes::ValueType::reference
2948    Data::getDataAtOffsetRW(DataTypes::ValueType::size_type i)
2949    {
2950        checkExclusiveWrite();
2951        return getReady()->getDataAtOffsetRW(i);
2952    }
2953    
2954    // This method is not thread-safe
2955  DataTypes::ValueType::const_reference  DataTypes::ValueType::const_reference
2956  Data::getDataAtOffset(DataTypes::ValueType::size_type i) const  Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i)
2957  {  {
2958      if (isLazy())      forceResolve();
2959      {      return getReady()->getDataAtOffsetRO(i);
     throw DataException("Programmer error - getDataAtOffset not permitted on lazy data (object is const which prevents resolving).");  
     }  
     return getReady()->getDataAtOffset(i);  
2960  }  }
2961    
2962    
2963  DataTypes::ValueType::reference  // DataTypes::ValueType::const_reference
2964  Data::getDataAtOffset(DataTypes::ValueType::size_type i)  // Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i) const
2965  {  // {
2966  //     if (isLazy())  //     if (isLazy())
2967  //     {  //     {
2968  //  throw DataException("getDataAtOffset not permitted on lazy data.");  //  throw DataException("Programmer error - getDataAtOffsetRO() not permitted on Lazy Data (object is const which prevents resolving).");
2969  //     }  //     }
2970      FORCERESOLVE;  //     return getReady()->getDataAtOffsetRO(i);
2971      return getReady()->getDataAtOffset(i);  // }
2972  }  
2973    
2974  DataTypes::ValueType::const_reference  DataTypes::ValueType::const_reference
2975  Data::getDataPoint(int sampleNo, int dataPointNo) const  Data::getDataPointRO(int sampleNo, int dataPointNo)
2976  {  {
2977      forceResolve();
2978    if (!isReady())    if (!isReady())
2979    {    {
2980      throw DataException("Programmer error - getDataPoint() not permitted on Lazy Data (object is const which prevents resolving).");      throw DataException("Programmer error -getDataPointRO() not permitted on Lazy Data.");
2981    }    }
2982    else    else
2983    {    {
2984      const DataReady* dr=getReady();      const DataReady* dr=getReady();
2985      return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));      return dr->getDataAtOffsetRO(dr->getPointOffset(sampleNo, dataPointNo));
2986    }    }
2987  }  }
2988    
2989    
2990    
2991    
2992  DataTypes::ValueType::reference  DataTypes::ValueType::reference
2993  Data::getDataPoint(int sampleNo, int dataPointNo)  Data::getDataPointRW(int sampleNo, int dataPointNo)
2994  {  {
2995    FORCERESOLVE;    checkExclusiveWrite();
2996    if (!isReady())    DataReady* dr=getReady();
2997    {    return dr->getDataAtOffsetRW(dr->getPointOffset(sampleNo, dataPointNo));
2998      throw DataException("Programmer error - getDataPoint() not permitted on Lazy Data.");  }
2999    }  
3000    else  BufferGroup*
3001    {  Data::allocSampleBuffer() const
3002      DataReady* dr=getReady();  {
3003      return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));       if (isLazy())
3004    }       {
3005        #ifdef _OPENMP
3006        int tnum=omp_get_max_threads();
3007        #else
3008        int tnum=1;
3009        #endif
3010        return new BufferGroup(getSampleBufferSize(),tnum);
3011         }
3012         else
3013         {
3014        return NULL;
3015         }
3016    }
3017    
3018    void
3019    Data::freeSampleBuffer(BufferGroup* bufferg)
3020    {
3021         if (bufferg!=0)
3022         {
3023        delete bufferg;
3024         }
3025    }
3026    
3027    
3028    Data
3029    Data::interpolateFromTable2DP(boost::python::object table, double Amin, double Astep,
3030            Data& B, double Bmin, double Bstep, double undef, bool check_boundaries)
3031    {
3032        WrappedArray t(table);
3033        return interpolateFromTable2D(t, Amin, Astep, undef, B, Bmin, Bstep,check_boundaries);
3034    }
3035    
3036    Data
3037    Data::interpolateFromTable1DP(boost::python::object table, double Amin, double Astep,
3038                              double undef,bool check_boundaries)
3039    {
3040        WrappedArray t(table);
3041        return interpolateFromTable1D(t, Amin, Astep, undef, check_boundaries);
3042    }
3043    
3044    
3045    Data
3046    Data::interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,
3047                     double undef, bool check_boundaries)
3048    {
3049        int error=0;
3050        if ((getDataPointRank()!=0))
3051        {
3052            throw DataException("Input to 1D interpolation must be scalar");
3053        }
3054        if (table.getRank()!=1)
3055        {
3056            throw DataException("Table for 1D interpolation must be 1D");
3057        }
3058        if (Astep<=0)
3059        {
3060            throw DataException("Astep must be positive");
3061        }
3062        if (!isExpanded())
3063        {
3064            expand();
3065        }
3066        Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);
3067        do                                   // to make breaks useful
3068        {
3069            try
3070            {
3071                int numpts=getNumDataPoints();
3072                const DataVector& adat=getReady()->getVectorRO();
3073                DataVector& rdat=res.getReady()->getVectorRW();
3074                int twidth=table.getShape()[0]-1;
3075                for (int l=0; l<numpts; ++l)
3076                {
3077                    double a=adat[l];
3078                    int x=static_cast<int>(((a-Amin)/Astep));
3079                                    if (check_boundaries) {
3080                        if ((a<Amin) || (x<0))
3081                        {
3082                            error=1;
3083                            break;  
3084                        }
3085                        if (a>Amin+Astep*twidth)
3086                        {
3087                            error=4;
3088                            break;
3089                        }
3090                                    }
3091                                    if (x<0) x=0;
3092                                    if (x>twidth) x=twidth;
3093    
3094                    if (x==twidth)          // value is on the far end of the table
3095                    {
3096                        double e=table.getElt(x);
3097                        if (e>undef)
3098                        {
3099                            error=2;
3100                            break;
3101                        }
3102                        rdat[l]=e;
3103                    }
3104                    else        // x and y are in bounds
3105                    {
3106                        double e=table.getElt(x);
3107                        double w=table.getElt(x+1);
3108                        if ((e>undef) || (w>undef))
3109                        {
3110                            error=2;
3111                            break;
3112                        }
3113                // map x*Astep <= a << (x+1)*Astep to [-1,1]
3114                        double la = 2.0*(a-Amin-(x*Astep))/Astep-1;
3115                        rdat[l]=((1-la)*e + (1+la)*w)/2;
3116                    }
3117                }
3118            } catch (DataException d)
3119            {
3120                error=3;
3121                break;
3122            }
3123        } while (false);
3124    #ifdef PASO_MPI
3125        int rerror=0;
3126        MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );
3127        error=rerror;
3128    #endif
3129        if (error)
3130        {
3131            switch (error)
3132            {
3133                case 1: throw DataException("Value below lower table range.");
3134                case 2: throw DataException("Interpolated value too large");
3135                case 4: throw DataException("Value greater than upper table range.");
3136                default:
3137                    throw DataException("Unknown error in interpolation");      
3138            }
3139        }
3140        return res;
3141    }
3142    
3143            
3144    Data
3145    Data::interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
3146                           double undef, Data& B, double Bmin, double Bstep, bool check_boundaries)
3147    {
3148        int error=0;
3149        if ((getDataPointRank()!=0) || (B.getDataPointRank()!=0))
3150        {
3151            throw DataException("Inputs to 2D interpolation must be scalar");
3152        }
3153        if (table.getRank()!=2)
3154        {
3155        throw DataException("Table for 2D interpolation must be 2D");
3156        }
3157        if ((Astep<=0) || (Bstep<=0))
3158        {
3159        throw DataException("Astep and Bstep must be postive");
3160        }
3161        if (getFunctionSpace()!=B.getFunctionSpace())
3162        {
3163        Data n=B.interpolate(getFunctionSpace());
3164        return interpolateFromTable2D(table, Amin, Astep, undef,
3165            n , Bmin, Bstep, check_boundaries);
3166        }
3167        if (!isExpanded())
3168        {
3169        expand();
3170        }
3171        if (!B.isExpanded())
3172        {
3173        B.expand();
3174        }
3175        Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);
3176        do                                   // to make breaks useful
3177        {
3178        try
3179        {
3180            int numpts=getNumDataPoints();
3181            const DataVector& adat=getReady()->getVectorRO();
3182            const DataVector& bdat=B.getReady()->getVectorRO();
3183            DataVector& rdat=res.getReady()->getVectorRW();
3184            const DataTypes::ShapeType& ts=table.getShape();
3185            int twx=ts[0]-1;    // table width x
3186            int twy=ts[1]-1;    // table width y
3187            for (int l=0; l<numpts; ++l)
3188            {
3189            double a=adat[l];
3190            double b=bdat[l];
3191            int x=static_cast<int>(((a-Amin)/Astep));
3192            int y=static_cast<int>(((b-Bmin)/Bstep));
3193                    if (check_boundaries) {
3194                    if ( (a<Amin) || (b<Bmin) || (x<0) || (y<0) )
3195                    {
3196                        error=1;
3197                        break;  
3198                    }
3199                    if ( (a>Amin+Astep*twx) || (b>Bmin+Bstep*twy) )
3200                    {
3201                        error=4;
3202                        break;
3203                    }
3204                    }
3205                    if (x<0) x=0;
3206                    if (y<0) y=0;
3207                    if (x>twx) x=twx;
3208                    if (y>twx) y=twy;
3209    
3210                    if (x == twx ) {
3211                         if (y == twy ) {
3212                     double sw=table.getElt(x,y);
3213                     if ((sw>undef))
3214                     {
3215                     error=2;
3216                     break;
3217                     }
3218                     rdat[l]=sw;
3219    
3220                         } else {
3221                     double sw=table.getElt(x,y);
3222                     double nw=table.getElt(x,y+1);
3223                     if ((sw>undef) || (nw>undef))
3224                     {
3225                     error=2;
3226                     break;
3227                     }
3228                     double lb = 2.0*(b-Bmin-(y*Bstep))/Bstep-1;
3229                     rdat[l]=((1-lb)*sw + (1+lb)*nw )/2.;
3230    
3231                         }
3232                    } else {
3233                         if (y == twy ) {
3234                     double sw=table.getElt(x,y);
3235                     double se=table.getElt(x+1,y);
3236                     if ((sw>undef) || (se>undef) )
3237                     {
3238                     error=2;
3239                     break;
3240                     }
3241                     double la = 2.0*(a-Amin-(x*Astep))/Astep-1;
3242                     rdat[l]=((1-la)*sw + (1+la)*se )/2;
3243    
3244                         } else {
3245                     double sw=table.getElt(x,y);
3246                     double nw=table.getElt(x,y+1);
3247                     double se=table.getElt(x+1,y);
3248                     double ne=table.getElt(x+1,y+1);
3249                     if ((sw>undef) || (nw>undef) || (se>undef) || (ne>undef))
3250                     {
3251                     error=2;
3252                     break;
3253                     }
3254                     // map x*Astep <= a << (x+1)*Astep to [-1,1]
3255                     // same with b
3256                     double la = 2.0*(a-Amin-(x*Astep))/Astep-1;
3257                     double lb = 2.0*(b-Bmin-(y*Bstep))/Bstep-1;
3258                     rdat[l]=((1-la)*(1-lb)*sw + (1-la)*(1+lb)*nw +
3259                          (1+la)*(1-lb)*se + (1+la)*(1+lb)*ne)/4;
3260                         }
3261                    }
3262            }
3263        } catch (DataException d)
3264        {
3265            error=3;
3266            break;
3267        }
3268        } while (false);
3269    #ifdef PASO_MPI
3270        int rerror=0;
3271        MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );
3272        error=rerror;
3273    #endif
3274        if (error)
3275        {
3276        switch (error)
3277        {
3278                case 1: throw DataException("Value below lower table range.");
3279                case 2: throw DataException("Interpolated value too large");
3280                case 4: throw DataException("Value greater than upper table range.");
3281                        default:
3282                          throw DataException("Unknown error in interpolation");        
3283        }
3284        }
3285        return res;
3286  }  }
3287    
3288    
# Line 3000  Data::print() Line 3298  Data::print()
3298    {    {
3299      printf( "[%6d]", i );      printf( "[%6d]", i );
3300      for( j=0; j<getNumDataPointsPerSample(); j++ )      for( j=0; j<getNumDataPointsPerSample(); j++ )
3301        printf( "\t%10.7g", (getSampleData(i))[j] );        printf( "\t%10.7g", (getSampleDataRW(i))[j] );    // doesn't really need RW access
3302      printf( "\n" );      printf( "\n" );
3303    }    }
3304  }  }
# Line 3020  Data::dump(const std::string fileName) c Line 3318  Data::dump(const std::string fileName) c
3318            return m_data->dump(fileName);            return m_data->dump(fileName);
3319      }      }
3320    }    }
3321    catch (exception& e)    catch (std::exception& e)
3322    {    {
3323          cout << e.what() << endl;          cout << e.what() << endl;
3324    }    }
# Line 3062  Data::get_MPIComm() const Line 3360  Data::get_MPIComm() const
3360  #endif  #endif
3361  }  }
3362    
   

Legend:
Removed from v.2037  
changed lines
  Added in v.2673

  ViewVC Help
Powered by ViewVC 1.1.26