/[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 2049 by phornby, Mon Nov 17 08:54:33 2008 UTC revision 2735 by jfenwick, Mon Nov 2 02:03:24 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    #define CHECK_DO_CRES escriptParams.getRESOLVE_COLLECTIVE()
91    
92    namespace
93    {
94    
95    template <class ARR>
96    inline
97    boost::python::tuple
98    pointToTuple1(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
99    {
100        using namespace boost::python;
101        using boost::python::tuple;
102        using boost::python::list;
103    
104        list l;
105        unsigned int dim=shape[0];
106        for (size_t i=0;i<dim;++i)
107        {
108        l.append(v[i+offset]);
109        }
110        return tuple(l);
111    }
112    
113    template <class ARR>
114    inline
115    boost::python::tuple
116    pointToTuple2(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
117    {
118        using namespace boost::python;
119        using boost::python::tuple;
120        using boost::python::list;
121    
122        unsigned int shape0=shape[0];
123        unsigned int shape1=shape[1];
124        list lj;
125        for (size_t j=0;j<shape0;++j)
126        {
127            list li;
128        for (size_t i=0;i<shape1;++i)
129        {
130            li.append(v[offset+DataTypes::getRelIndex(shape,j,i)]);
131        }
132        lj.append(tuple(li));
133        }
134        return tuple(lj);
135    }
136    
137    template <class ARR>
138    inline
139    boost::python::tuple
140    pointToTuple3(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
141    {
142        using namespace boost::python;
143        using boost::python::tuple;
144        using boost::python::list;
145    
146        unsigned int shape0=shape[0];
147        unsigned int shape1=shape[1];
148        unsigned int shape2=shape[2];
149    
150        list lk;
151        for (size_t k=0;k<shape0;++k)
152        {
153            list lj;
154        for (size_t j=0;j<shape1;++j)
155        {
156            list li;
157            for (size_t i=0;i<shape2;++i)
158            {
159                    li.append(v[offset+DataTypes::getRelIndex(shape,k,j,i)]);
160                }
161            lj.append(tuple(li));
162            }
163            lk.append(tuple(lj));
164        }
165        return tuple(lk);
166    }
167    
168    template <class ARR>
169    inline
170    boost::python::tuple
171    pointToTuple4(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
172    {
173        using namespace boost::python;
174        using boost::python::tuple;
175        using boost::python::list;
176    
177        unsigned int shape0=shape[0];
178        unsigned int shape1=shape[1];
179        unsigned int shape2=shape[2];
180        unsigned int shape3=shape[3];
181    
182        list ll;
183        for (size_t l=0;l<shape0;++l)
184        {
185            list lk;
186        for (size_t k=0;k<shape1;++k)
187        {
188                list lj;
189                for (size_t j=0;j<shape2;++j)
190                {
191                    list li;
192                    for (size_t i=0;i<shape3;++i)
193                    {
194                        li.append(v[offset+DataTypes::getRelIndex(shape,l,k,j,i)]);
195                    }
196                    lj.append(tuple(li));
197                }
198                lk.append(tuple(lj));
199        }
200            ll.append(tuple(lk));
201        }
202        return tuple(ll);
203    }
204    
205    
206    // This should be safer once the DataC RO changes have been brought in
207    template <class ARR>
208    boost::python::tuple
209    pointToTuple( const DataTypes::ShapeType& shape,ARR v)
210    {
211       using namespace boost::python;
212       using boost::python::list;
213       int rank=shape.size();
214       if (rank==0)
215       {
216        return make_tuple(v[0]);
217       }
218       else if (rank==1)
219       {
220            return pointToTuple1(shape,v,0);
221       }
222       else if (rank==2)
223       {
224        return pointToTuple2(shape,v,0);
225       }
226       else if (rank==3)
227       {
228        return pointToTuple3(shape,v,0);
229       }
230       else if (rank==4)
231       {
232        return pointToTuple4(shape,v,0);
233       }
234       else
235         throw DataException("Unknown rank in pointToTuple.");
236    }
237    
238    }  // anonymous namespace
239    
240  Data::Data()  Data::Data()
241        : m_shared(false), m_lazy(false)
242  {  {
243    //    //
244    // Default data is type DataEmpty    // Default data is type DataEmpty
245    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
246    m_data=temp->getPtr();  //   m_data=temp->getPtr();
247      set_m_data(temp->getPtr());
248    m_protected=false;    m_protected=false;
249  }  }
250    
# Line 60  Data::Data(double value, Line 252  Data::Data(double value,
252             const tuple& shape,             const tuple& shape,
253             const FunctionSpace& what,             const FunctionSpace& what,
254             bool expanded)             bool expanded)
255        : m_shared(false), m_lazy(false)
256  {  {
257    DataTypes::ShapeType dataPointShape;    DataTypes::ShapeType dataPointShape;
258    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 269  Data::Data(double value,
269         const DataTypes::ShapeType& dataPointShape,         const DataTypes::ShapeType& dataPointShape,
270         const FunctionSpace& what,         const FunctionSpace& what,
271             bool expanded)             bool expanded)
272        : m_shared(false), m_lazy(false)
273  {  {
274    int len = DataTypes::noValues(dataPointShape);    int len = DataTypes::noValues(dataPointShape);
   
275    DataVector temp_data(len,value,len);    DataVector temp_data(len,value,len);
 //   DataArrayView temp_dataView(temp_data, dataPointShape);  
   
 //   initialise(temp_dataView, what, expanded);  
276    initialise(temp_data, dataPointShape, what, expanded);    initialise(temp_data, dataPointShape, what, expanded);
   
277    m_protected=false;    m_protected=false;
278  }  }
279    
280  Data::Data(const Data& inData)  Data::Data(const Data& inData)
281        : m_shared(false), m_lazy(false)
282  {  {
283    m_data=inData.m_data;    set_m_data(inData.m_data);
284    m_protected=inData.isProtected();    m_protected=inData.isProtected();
285  }  }
286    
287    
288  Data::Data(const Data& inData,  Data::Data(const Data& inData,
289             const DataTypes::RegionType& region)             const DataTypes::RegionType& region)
290        : m_shared(false), m_lazy(false)
291  {  {
292    DataAbstract_ptr dat=inData.m_data;    DataAbstract_ptr dat=inData.m_data;
293    if (inData.isLazy())    if (inData.isLazy())
# Line 110  Data::Data(const Data& inData, Line 301  Data::Data(const Data& inData,
301    //    //
302    // Create Data which is a slice of another Data    // Create Data which is a slice of another Data
303    DataAbstract* tmp = dat->getSlice(region);    DataAbstract* tmp = dat->getSlice(region);
304    m_data=DataAbstract_ptr(tmp);    set_m_data(DataAbstract_ptr(tmp));
305    m_protected=false;    m_protected=false;
306    
307  }  }
308    
309  Data::Data(const Data& inData,  Data::Data(const Data& inData,
310             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
311        : m_shared(false), m_lazy(false)
312  {  {
313    if (inData.isEmpty())    if (inData.isEmpty())
314    {    {
315      throw DataException("Error - will not interpolate for instances of DataEmpty.");      throw DataException("Error - will not interpolate for instances of DataEmpty.");
316    }    }
317    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
318      m_data=inData.m_data;      set_m_data(inData.m_data);
319    }    }
320    else    else
321    {    {
# Line 130  Data::Data(const Data& inData, Line 323  Data::Data(const Data& inData,
323      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
324        if (!inData.probeInterpolation(functionspace))        if (!inData.probeInterpolation(functionspace))
325        {           // 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
326      throw FunctionSpaceException("Call to probeInterpolation returned false for DataConstant.");      throw FunctionSpaceException("Cannot interpolate across to the domain of the specified FunctionSpace. (DataConstant)");
327        }        }
328        // 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
329        DataReady_ptr dr=inData.m_data->resolve();        DataReady_ptr dr=inData.m_data->resolve();
330        DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVector());          DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVectorRO());
331        m_data=DataAbstract_ptr(dc);  //       m_data=DataAbstract_ptr(dc);
332          set_m_data(DataAbstract_ptr(dc));
333      } else {      } else {
334        Data tmp(0,inData.getDataPointShape(),functionspace,true);        Data tmp(0,inData.getDataPointShape(),functionspace,true);
335        // 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 343  Data::Data(const Data& inData,
343        } else {        } else {
344          inDataDomain->interpolateACross(tmp,inData);          inDataDomain->interpolateACross(tmp,inData);
345        }        }
346        m_data=tmp.m_data;  //       m_data=tmp.m_data;
347          set_m_data(tmp.m_data);
348      }      }
349    }    }
350    m_protected=false;    m_protected=false;
351  }  }
352    
353  Data::Data(DataAbstract* underlyingdata)  Data::Data(DataAbstract* underlyingdata)
354        : m_shared(false), m_lazy(false)
355  {  {
356  //  m_data=shared_ptr<DataAbstract>(underlyingdata);      set_m_data(underlyingdata->getPtr());
     m_data=underlyingdata->getPtr();  
357      m_protected=false;      m_protected=false;
358  }  }
359    
360  Data::Data(DataAbstract_ptr underlyingdata)  Data::Data(DataAbstract_ptr underlyingdata)
361        : m_shared(false), m_lazy(false)
362  {  {
363      m_data=underlyingdata;      set_m_data(underlyingdata);
364      m_protected=false;      m_protected=false;
365  }  }
366    
   
 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;  
 }*/  
   
367  Data::Data(const DataTypes::ValueType& value,  Data::Data(const DataTypes::ValueType& value,
368           const DataTypes::ShapeType& shape,           const DataTypes::ShapeType& shape,
369                   const FunctionSpace& what,                   const FunctionSpace& what,
370                   bool expanded)                   bool expanded)
371        : m_shared(false), m_lazy(false)
372  {  {
373     initialise(value,shape,what,expanded);     initialise(value,shape,what,expanded);
374     m_protected=false;     m_protected=false;
# Line 198  Data::Data(const DataTypes::ValueType& v Line 378  Data::Data(const DataTypes::ValueType& v
378  Data::Data(const object& value,  Data::Data(const object& value,
379         const FunctionSpace& what,         const FunctionSpace& what,
380             bool expanded)             bool expanded)
381        : m_shared(false), m_lazy(false)
382  {  {
383    numeric::array asNumArray(value);    WrappedArray w(value);
384    initialise(asNumArray,what,expanded);    initialise(w,what,expanded);
385    m_protected=false;    m_protected=false;
386  }  }
387    
388    
389  Data::Data(const object& value,  Data::Data(const object& value,
390             const Data& other)             const Data& other)
391        : m_shared(false), m_lazy(false)
392  {  {
393    numeric::array asNumArray(value);    WrappedArray w(value);
   
   // extract the shape of the numarray  
   DataTypes::ShapeType tempShape=DataTypes::shapeFromNumArray(asNumArray);  
 // /*  for (int i=0; i < asNumArray.getrank(); i++) {  
 //     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);  
394    
395    //    // extract the shape of the array
396    // Create DataConstant using the given value and all other parameters    const DataTypes::ShapeType& tempShape=w.getShape();
397    // copied from other. If value is a rank 0 object this Data    if (w.getRank()==0) {
   // will assume the point data shape of other.  
   
   if (DataTypes::getRank(tempShape)/*temp_dataView.getRank()*/==0) {  
398    
399    
400      // get the space for the data vector      // get the space for the data vector
401      int len1 = DataTypes::noValues(tempShape);      int len1 = DataTypes::noValues(tempShape);
402      DataVector temp_data(len1, 0.0, len1);      DataVector temp_data(len1, 0.0, len1);
403      temp_data.copyFromNumArray(asNumArray);      temp_data.copyFromArray(w,1);
404    
405      int len = DataTypes::noValues(other.getDataPointShape());      int len = DataTypes::noValues(other.getDataPointShape());
406    
407      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);  
   
408      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
409  //     boost::shared_ptr<DataAbstract> sp(t);  //     m_data=DataAbstract_ptr(t);
410  //     m_data=sp;      set_m_data(DataAbstract_ptr(t));
     m_data=DataAbstract_ptr(t);  
411    
412    } else {    } else {
413      //      //
414      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
415  //     initialise(temp_dataView, other.getFunctionSpace(), false);      DataConstant* t=new DataConstant(w,other.getFunctionSpace());
416      DataConstant* t=new DataConstant(asNumArray,other.getFunctionSpace());  //     m_data=DataAbstract_ptr(t);
417  //     boost::shared_ptr<DataAbstract> sp(t);      set_m_data(DataAbstract_ptr(t));
 //     m_data=sp;  
     m_data=DataAbstract_ptr(t);  
418    }    }
419    m_protected=false;    m_protected=false;
420  }  }
421    
422  Data::~Data()  Data::~Data()
423  {  {
424      set_m_data(DataAbstract_ptr());
425  }  }
426    
427    
428    // only call in thread safe contexts.
429    // This method should be atomic
430    void Data::set_m_data(DataAbstract_ptr p)
431    {
432      if (m_data.get()!=0)  // release old ownership
433      {
434        m_data->removeOwner(this);
435      }
436      if (p.get()!=0)
437      {
438        m_data=p;
439        m_data->addOwner(this);
440        m_shared=m_data->isShared();
441        m_lazy=m_data->isLazy();
442      }
443    }
444    
445  void  void Data::initialise(const WrappedArray& value,
 Data::initialise(const boost::python::numeric::array& value,  
446                   const FunctionSpace& what,                   const FunctionSpace& what,
447                   bool expanded)                   bool expanded)
448  {  {
# Line 277  Data::initialise(const boost::python::nu Line 453  Data::initialise(const boost::python::nu
453    // within the shared_ptr constructor.    // within the shared_ptr constructor.
454    if (expanded) {    if (expanded) {
455      DataAbstract* temp=new DataExpanded(value, what);      DataAbstract* temp=new DataExpanded(value, what);
456  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
457  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
458    } else {    } else {
459      DataAbstract* temp=new DataConstant(value, what);      DataAbstract* temp=new DataConstant(value, what);
460  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
461  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
462    }    }
463  }  }
464    
# Line 302  Data::initialise(const DataTypes::ValueT Line 476  Data::initialise(const DataTypes::ValueT
476    // within the shared_ptr constructor.    // within the shared_ptr constructor.
477    if (expanded) {    if (expanded) {
478      DataAbstract* temp=new DataExpanded(what, shape, value);      DataAbstract* temp=new DataExpanded(what, shape, value);
479  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
480  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
481    } else {    } else {
482      DataAbstract* temp=new DataConstant(what, shape, value);      DataAbstract* temp=new DataConstant(what, shape, value);
483  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
484  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
485    }    }
486  }  }
487    
488    
 // 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";  
 //  }  
 // }  
   
489  escriptDataC  escriptDataC
490  Data::getDataC()  Data::getDataC()
491  {  {
# Line 410  Data::getDataC() const Line 502  Data::getDataC() const
502    return temp;    return temp;
503  }  }
504    
505    size_t
506    Data::getSampleBufferSize() const
507    {
508      return m_data->getSampleBufferSize();
509    }
510    
511    
512  const boost::python::tuple  const boost::python::tuple
513  Data::getShapeTuple() const  Data::getShapeTuple() const
514  {  {
# Line 435  Data::getShapeTuple() const Line 534  Data::getShapeTuple() const
534  // 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.
535  // 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
536  // 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
537  Data*  Data
538  Data::copySelf()  Data::copySelf()
539  {  {
540     DataAbstract* temp=m_data->deepCopy();     DataAbstract* temp=m_data->deepCopy();
541     return new Data(temp);     return Data(temp);
542  }  }
543    
544  void  void
# Line 447  Data::copy(const Data& other) Line 546  Data::copy(const Data& other)
546  {  {
547    DataAbstract* temp=other.m_data->deepCopy();    DataAbstract* temp=other.m_data->deepCopy();
548    DataAbstract_ptr p=temp->getPtr();    DataAbstract_ptr p=temp->getPtr();
549    m_data=p;  //   m_data=p;
550      set_m_data(p);
551  }  }
552    
553    
554  Data  Data
555  Data::delay()  Data::delay()
556  {  {
557    DataLazy* dl=new DataLazy(m_data);    if (!isLazy())
558    return Data(dl);    {
559          DataLazy* dl=new DataLazy(m_data);
560          return Data(dl);
561      }
562      return *this;
563  }  }
564    
565  void  void
# Line 463  Data::delaySelf() Line 567  Data::delaySelf()
567  {  {
568    if (!isLazy())    if (!isLazy())
569    {    {
570      m_data=(new DataLazy(m_data))->getPtr();  //  m_data=(new DataLazy(m_data))->getPtr();
571        set_m_data((new DataLazy(m_data))->getPtr());
572    }    }
573  }  }
574    
575    
576    // For lazy data, it would seem that DataTagged will need to be treated differently since even after setting all tags
577    // to zero, all the tags from all the DataTags would be in the result.
578    // However since they all have the same value (0) whether they are there or not should not matter.
579    // So I have decided that for all types this method will create a constant 0.
580    // It can be promoted up as required.
581    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
582    // but we can deal with that if it arises.
583    //
584  void  void
585  Data::setToZero()  Data::setToZero()
586  {  {
# Line 474  Data::setToZero() Line 588  Data::setToZero()
588    {    {
589       throw DataException("Error - Operations not permitted on instances of DataEmpty.");       throw DataException("Error - Operations not permitted on instances of DataEmpty.");
590    }    }
591    m_data->setToZero();    if (isLazy())
592      {
593         DataTypes::ValueType v(getNoValues(),0);
594         DataConstant* dc=new DataConstant(getFunctionSpace(),getDataPointShape(),v);
595         DataLazy* dl=new DataLazy(dc->getPtr());
596         set_m_data(dl->getPtr());
597      }
598      else
599      {
600         exclusiveWrite();
601         m_data->setToZero();
602      }
603  }  }
604    
605    
606  void  void
607  Data::copyWithMask(const Data& other,  Data::copyWithMask(const Data& other,
608                     const Data& mask)                     const Data& mask)
# Line 538  Data::copyWithMask(const Data& other, Line 664  Data::copyWithMask(const Data& other,
664    {    {
665      throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");      throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
666    }    }
667      unsigned int selfrank=getDataPointRank();
668      unsigned int otherrank=other2.getDataPointRank();
669      unsigned int maskrank=mask2.getDataPointRank();
670      if ((selfrank==0) && (otherrank>0 || maskrank>0))
671      {
672        // to get here we must be copying from a large object into a scalar
673        // I am not allowing this.
674        // If you are calling copyWithMask then you are considering keeping some existing values
675        // and so I'm going to assume that you don't want your data objects getting a new shape.
676        throw DataException("Attempt to copyWithMask into a scalar from an object or mask with rank>0.");
677      }
678      exclusiveWrite();
679    // Now we iterate over the elements    // Now we iterate over the elements
680    DataVector& self=getReadyPtr()->getVector();    DataVector& self=getReady()->getVectorRW();;
681    const DataVector& ovec=other2.getReadyPtr()->getVector();    const DataVector& ovec=other2.getReadyPtr()->getVectorRO();
682    const DataVector& mvec=mask2.getReadyPtr()->getVector();    const DataVector& mvec=mask2.getReadyPtr()->getVectorRO();
683    if ((self.size()!=ovec.size()) || (self.size()!=mvec.size()))  
684      if ((selfrank>0) && (otherrank==0) &&(maskrank==0))
685      {
686        // Not allowing this combination.
687        // it is not clear what the rank of the target should be.
688        // Should it be filled with the scalar (rank stays the same);
689        // or should the target object be reshaped to be a scalar as well.
690        throw DataException("Attempt to copyWithMask from scalar mask and data into non-scalar target.");
691      }
692      if ((selfrank>0) && (otherrank>0) &&(maskrank==0))
693      {
694        if (mvec[0]>0)      // copy whole object if scalar is >0
695        {
696            copy(other);
697        }
698        return;
699      }
700      if (isTagged())       // so all objects involved will also be tagged
701    {    {
702      throw DataException("Error - size mismatch in arguments to copyWithMask.");      // note the !
703        if (!((getDataPointShape()==mask2.getDataPointShape()) &&
704            ((other2.getDataPointShape()==mask2.getDataPointShape()) || (otherrank==0))))
705        {
706            throw DataException("copyWithMask, shape mismatch.");
707        }
708    
709        // We need to consider the possibility that tags are missing or in the wrong order
710        // My guiding assumption here is: All tagged Datas are assumed to have the default value for
711        // all tags which are not explicitly defined
712    
713        const DataTagged* mptr=dynamic_cast<const DataTagged*>(mask2.m_data.get());
714        const DataTagged* optr=dynamic_cast<const DataTagged*>(other2.m_data.get());
715        DataTagged* tptr=dynamic_cast<DataTagged*>(m_data.get());
716    
717        // first, add any tags missing from other or mask
718        const DataTagged::DataMapType& olookup=optr->getTagLookup();
719            const DataTagged::DataMapType& mlookup=mptr->getTagLookup();
720        const DataTagged::DataMapType& tlookup=tptr->getTagLookup();
721        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
722        for (i=olookup.begin();i!=olookup.end();i++)
723        {
724               tptr->addTag(i->first);
725            }
726            for (i=mlookup.begin();i!=mlookup.end();i++) {
727               tptr->addTag(i->first);
728            }
729        // now we know that *this has all the required tags but they aren't guaranteed to be in
730        // the same order
731    
732        // There are two possibilities: 1. all objects have the same rank. 2. other is a scalar
733        if ((selfrank==otherrank) && (otherrank==maskrank))
734        {
735            for (i=tlookup.begin();i!=tlookup.end();i++)
736            {
737                // get the target offset
738                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
739                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
740                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
741                for (int j=0;j<getDataPointSize();++j)
742                {
743                    if (mvec[j+moff]>0)
744                    {
745                        self[j+toff]=ovec[j+ooff];
746                    }
747                }
748                }
749            // now for the default value
750            for (int j=0;j<getDataPointSize();++j)
751            {
752                if (mvec[j+mptr->getDefaultOffset()]>0)
753                {
754                    self[j+tptr->getDefaultOffset()]=ovec[j+optr->getDefaultOffset()];
755                }
756            }
757        }
758        else    // other is a scalar
759        {
760            for (i=tlookup.begin();i!=tlookup.end();i++)
761            {
762                // get the target offset
763                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
764                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
765                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
766                for (int j=0;j<getDataPointSize();++j)
767                {
768                    if (mvec[j+moff]>0)
769                    {
770                        self[j+toff]=ovec[ooff];
771                    }
772                }
773                }
774            // now for the default value
775            for (int j=0;j<getDataPointSize();++j)
776            {
777                if (mvec[j+mptr->getDefaultOffset()]>0)
778                {
779                    self[j+tptr->getDefaultOffset()]=ovec[0];
780                }
781            }
782        }
783    
784        return;         // ugly
785      }
786      // mixed scalar and non-scalar operation
787      if ((selfrank>0) && (otherrank==0) && (mask2.getDataPointShape()==getDataPointShape()))
788      {
789            size_t num_points=self.size();
790        // OPENMP 3.0 allows unsigned loop vars.
791        #if defined(_OPENMP) && (_OPENMP < 200805)
792        long i;
793        #else
794        size_t i;
795        #endif  
796        size_t psize=getDataPointSize();    
797        #pragma omp parallel for private(i) schedule(static)
798        for (i=0;i<num_points;++i)
799        {
800            if (mvec[i]>0)
801            {
802                self[i]=ovec[i/psize];      // since this is expanded there is one scalar
803            }                   // dest point
804        }
805        return;         // ugly!
806      }
807      // tagged data is already taken care of so we only need to worry about shapes
808      // special cases with scalars are already dealt with so all we need to worry about is shape
809      if ((getDataPointShape()!=other2.getDataPointShape()) || getDataPointShape()!=mask2.getDataPointShape())
810      {
811        ostringstream oss;
812        oss <<"Error - size mismatch in arguments to copyWithMask.";
813        oss << "\nself_shape=" << DataTypes::shapeToString(getDataPointShape());
814        oss << " other2_shape=" << DataTypes::shapeToString(other2.getDataPointShape());
815        oss << " mask2_shape=" << DataTypes::shapeToString(mask2.getDataPointShape());
816        throw DataException(oss.str());
817    }    }
818    size_t num_points=self.size();    size_t num_points=self.size();
819    
# Line 564  Data::copyWithMask(const Data& other, Line 833  Data::copyWithMask(const Data& other,
833    }    }
834  }  }
835    
   
   
836  bool  bool
837  Data::isExpanded() const  Data::isExpanded() const
838  {  {
# Line 574  Data::isExpanded() const Line 841  Data::isExpanded() const
841  }  }
842    
843  bool  bool
844    Data::actsExpanded() const
845    {
846      return m_data->actsExpanded();
847    
848    }
849    
850    
851    bool
852  Data::isTagged() const  Data::isTagged() const
853  {  {
854    DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());    DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
# Line 597  Data::isConstant() const Line 872  Data::isConstant() const
872  bool  bool
873  Data::isLazy() const  Data::isLazy() const
874  {  {
875    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
876  }  }
877    
878  // 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 903  Data::expand()
903    if (isConstant()) {    if (isConstant()) {
904      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
905      DataAbstract* temp=new DataExpanded(*tempDataConst);      DataAbstract* temp=new DataExpanded(*tempDataConst);
906  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
907  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
908    } else if (isTagged()) {    } else if (isTagged()) {
909      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
910      DataAbstract* temp=new DataExpanded(*tempDataTag);      DataAbstract* temp=new DataExpanded(*tempDataTag);
911  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
912  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
913    } else if (isExpanded()) {    } else if (isExpanded()) {
914      //      //
915      // do nothing      // do nothing
# Line 656  Data::tag() Line 929  Data::tag()
929    if (isConstant()) {    if (isConstant()) {
930      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
931      DataAbstract* temp=new DataTagged(*tempDataConst);      DataAbstract* temp=new DataTagged(*tempDataConst);
932  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
933  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
934    } else if (isTagged()) {    } else if (isTagged()) {
935      // do nothing      // do nothing
936    } else if (isExpanded()) {    } else if (isExpanded()) {
# Line 671  Data::tag() Line 943  Data::tag()
943       {       {
944      throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");      throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");
945       }       }
946       m_data=res;      //      m_data=res;
947         set_m_data(res);
948       tag();       tag();
949    } else {    } else {
950      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 956  Data::resolve()
956  {  {
957    if (isLazy())    if (isLazy())
958    {    {
959       m_data=m_data->resolve();  //      m_data=m_data->resolve();
960        set_m_data(m_data->resolve());
961    }    }
962  }  }
963    
964    void
965    Data::requireWrite()
966    {
967      resolve();
968      exclusiveWrite();
969    }
970    
971  Data  Data
972  Data::oneOver() const  Data::oneOver() const
973  {  {
974    if (isLazy())    MAKELAZYOP(RECIP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),RECIP);  
     return Data(c);  
   }  
975    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
976  }  }
977    
978  Data  Data
979  Data::wherePositive() const  Data::wherePositive() const
980  {  {
981    if (isLazy())    MAKELAZYOP(GZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GZ);  
     return Data(c);  
   }  
982    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
983  }  }
984    
985  Data  Data
986  Data::whereNegative() const  Data::whereNegative() const
987  {  {
988    if (isLazy())    MAKELAZYOP(LZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LZ);  
     return Data(c);  
   }  
989    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
990  }  }
991    
992  Data  Data
993  Data::whereNonNegative() const  Data::whereNonNegative() const
994  {  {
995    if (isLazy())    MAKELAZYOP(GEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GEZ);  
     return Data(c);  
   }  
996    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
997  }  }
998    
999  Data  Data
1000  Data::whereNonPositive() const  Data::whereNonPositive() const
1001  {  {
1002    if (isLazy())    MAKELAZYOP(LEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LEZ);  
     return Data(c);  
   }  
1003    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
1004  }  }
1005    
1006  Data  Data
1007  Data::whereZero(double tol) const  Data::whereZero(double tol) const
1008  {  {
1009    Data dataAbs=abs();  //   Data dataAbs=abs();
1010    return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
1011       MAKELAZYOPOFF(EZ,tol)
1012       return C_TensorUnaryOperation(*this, bind2nd(AbsLTE(),tol));
1013    
1014  }  }
1015    
1016  Data  Data
1017  Data::whereNonZero(double tol) const  Data::whereNonZero(double tol) const
1018  {  {
1019    Data dataAbs=abs();  //   Data dataAbs=abs();
1020    return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
1021      MAKELAZYOPOFF(NEZ,tol)
1022      return C_TensorUnaryOperation(*this, bind2nd(AbsGT(),tol));
1023    
1024  }  }
1025    
1026  Data  Data
# Line 767  bool Line 1033  bool
1033  Data::probeInterpolation(const FunctionSpace& functionspace) const  Data::probeInterpolation(const FunctionSpace& functionspace) const
1034  {  {
1035    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());  
 //     }  
 //   }  
1036  }  }
1037    
1038  Data  Data
# Line 813  Data::getDataPointSize() const Line 1069  Data::getDataPointSize() const
1069    return m_data->getNoValues();    return m_data->getNoValues();
1070  }  }
1071    
1072    
1073  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1074  Data::getLength() const  Data::getLength() const
1075  {  {
1076    return m_data->getLength();    return m_data->getLength();
1077  }  }
1078    
 const  
 boost::python::numeric::array  
 Data:: getValueOfDataPoint(int dataPointNo)  
 {  
   int i, j, k, l;  
   
   FORCERESOLVE;  
1079    
1080    //  // There is no parallelism here ... elements need to be added in the correct order.
1081    // determine the rank and shape of each data point  //   If we could presize the list and then fill in the elements it might work
1082    int dataPointRank = getDataPointRank();  //   This would need setting elements to be threadsafe.
1083    const DataTypes::ShapeType& dataPointShape = getDataPointShape();  //   Having mulitple C threads calling into one interpreter is aparently a no-no.
1084    const boost::python::object
1085    //  Data::toListOfTuples(bool scalarastuple)
1086    // create the numeric array to be returned  {
1087    boost::python::numeric::array numArray(0.0);      using namespace boost::python;
1088        using boost::python::list;
1089    //      if (get_MPISize()>1)
1090    // the shape of the returned numeric array will be the same      {
1091    // as that of the data point          throw DataException("::toListOfTuples is not available for MPI with more than one process.");
1092    int arrayRank = dataPointRank;      }
1093    const DataTypes::ShapeType& arrayShape = dataPointShape;      unsigned int rank=getDataPointRank();
1094        unsigned int size=getDataPointSize();
1095    
1096    //      int npoints=getNumDataPoints();
1097    // resize the numeric array to the shape just calculated      expand();           // This will also resolve if required
1098    if (arrayRank==0) {      const DataTypes::ValueType& vec=getReady()->getVectorRO();
1099      numArray.resize(1);      boost::python::list temp;
1100    }      temp.append(object());
1101    if (arrayRank==1) {      boost::python::list res(temp*npoints);// presize the list by the "[None] * npoints"  trick
1102      numArray.resize(arrayShape[0]);      if (rank==0)
1103    }      {
1104    if (arrayRank==2) {          long count;
1105      numArray.resize(arrayShape[0],arrayShape[1]);          if (scalarastuple)
1106    }          {
1107    if (arrayRank==3) {              for (count=0;count<npoints;++count)
1108      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);              {
1109    }          res[count]=make_tuple(vec[count]);
1110    if (arrayRank==4) {              }
1111      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);          }
1112    }          else
1113            {
1114                for (count=0;count<npoints;++count)
1115                {
1116                    res[count]=vec[count];
1117                }
1118            }
1119        }
1120        else if (rank==1)
1121        {
1122            size_t count;
1123            size_t offset=0;
1124            for (count=0;count<npoints;++count,offset+=size)
1125            {
1126                res[count]=pointToTuple1(getDataPointShape(), vec, offset);
1127            }
1128        }
1129        else if (rank==2)
1130        {
1131            size_t count;
1132            size_t offset=0;
1133            for (count=0;count<npoints;++count,offset+=size)
1134            {
1135            res[count]=pointToTuple2(getDataPointShape(), vec, offset);
1136            }
1137        }
1138        else if (rank==3)
1139        {
1140            size_t count;
1141            size_t offset=0;
1142            for (count=0;count<npoints;++count,offset+=size)
1143            {
1144                res[count]=pointToTuple3(getDataPointShape(), vec, offset);
1145            }
1146        }
1147        else if (rank==4)
1148        {
1149            size_t count;
1150            size_t offset=0;
1151            for (count=0;count<npoints;++count,offset+=size)
1152            {
1153                res[count]=pointToTuple4(getDataPointShape(), vec, offset);
1154            }
1155        }
1156        else
1157        {
1158            throw DataException("Unknown rank in ::toListOfTuples()");
1159        }
1160        return res;
1161    }
1162    
1163    if (getNumDataPointsPerSample()>0) {  const boost::python::object
1164    Data::getValueOfDataPointAsTuple(int dataPointNo)
1165    {
1166       forceResolve();
1167       if (getNumDataPointsPerSample()>0) {
1168         int sampleNo = dataPointNo/getNumDataPointsPerSample();         int sampleNo = dataPointNo/getNumDataPointsPerSample();
1169         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1170         //         //
1171         // Check a valid sample number has been supplied         // Check a valid sample number has been supplied
1172         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1173             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");             throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.");
1174         }         }
1175    
1176         //         //
1177         // Check a valid data point number has been supplied         // Check a valid data point number has been supplied
1178         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1179             throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");             throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid dataPointNoInSample.");
1180         }         }
1181         // TODO: global error handling         // TODO: global error handling
        // create a view of the data if it is stored locally  
 //       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);  
1182         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1183           return pointToTuple(getDataPointShape(),&(getDataAtOffsetRO(offset)));
1184      }
1185         switch( dataPointRank ){    else
1186              case 0 :    {
1187                  numArray[0] = getDataAtOffset(offset);      // The pre-numpy method would return an empty array of the given shape
1188                  break;      // I'm going to throw an exception because if we have zero points per sample we have problems
1189              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;  
     }  
1190    }    }
   //  
   // return the array  
   return numArray;  
1191    
1192  }  }
1193    
1194    
1195  void  void
1196  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
1197  {  {
1198      // this will throw if the value cannot be represented      // this will throw if the value cannot be represented
1199      boost::python::numeric::array num_array(py_object);      setValueOfDataPointToArray(dataPointNo,py_object);
     setValueOfDataPointToArray(dataPointNo,num_array);  
1200  }  }
1201    
1202  void  void
1203  Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)  Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::object& obj)
1204  {  {
1205    if (isProtected()) {    if (isProtected()) {
1206          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1207    }    }
1208    FORCERESOLVE;  
1209      WrappedArray w(obj);
1210    //    //
1211    // check rank    // check rank
1212    if (static_cast<unsigned int>(num_array.getrank())<getDataPointRank())    if (static_cast<unsigned int>(w.getRank())<getDataPointRank())
1213        throw DataException("Rank of numarray does not match Data object rank");        throw DataException("Rank of array does not match Data object rank");
1214    
1215    //    //
1216    // check shape of num_array    // check shape of array
1217    for (unsigned int i=0; i<getDataPointRank(); i++) {    for (unsigned int i=0; i<getDataPointRank(); i++) {
1218      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])      if (w.getShape()[i]!=getDataPointShape()[i])
1219         throw DataException("Shape of numarray does not match Data object rank");         throw DataException("Shape of array does not match Data object rank");
1220    }    }
1221    
1222      exclusiveWrite();
1223    
1224    //    //
1225    // make sure data is expanded:    // make sure data is expanded:
1226    //    //
# Line 949  Data::setValueOfDataPointToArray(int dat Line 1230  Data::setValueOfDataPointToArray(int dat
1230    if (getNumDataPointsPerSample()>0) {    if (getNumDataPointsPerSample()>0) {
1231         int sampleNo = dataPointNo/getNumDataPointsPerSample();         int sampleNo = dataPointNo/getNumDataPointsPerSample();
1232         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1233         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,w);
1234    } else {    } else {
1235         m_data->copyToDataPoint(-1, 0,num_array);         m_data->copyToDataPoint(-1, 0,w);
1236    }    }
1237  }  }
1238    
# Line 963  Data::setValueOfDataPoint(int dataPointN Line 1244  Data::setValueOfDataPoint(int dataPointN
1244    }    }
1245    //    //
1246    // make sure data is expanded:    // make sure data is expanded:
1247    FORCERESOLVE;    exclusiveWrite();
1248    if (!isExpanded()) {    if (!isExpanded()) {
1249      expand();      expand();
1250    }    }
# Line 977  Data::setValueOfDataPoint(int dataPointN Line 1258  Data::setValueOfDataPoint(int dataPointN
1258  }  }
1259    
1260  const  const
1261  boost::python::numeric::array  boost::python::object
1262  Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)  Data::getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo)
1263  {  {
1264    size_t length=0;    // This could be lazier than it is now
1265    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();  
1266    
1267    //    // copy datapoint into a buffer
1268    // create the numeric array to be returned    // broadcast buffer to all nodes
1269    boost::python::numeric::array numArray(0.0);    // convert buffer to tuple
1270      // return tuple
   //  
   // the shape of the returned numeric array will be the same  
   // as that of the data point  
   int arrayRank = dataPointRank;  
   const DataTypes::ShapeType& arrayShape = dataPointShape;  
1271    
1272    //    const DataTypes::ShapeType& dataPointShape = getDataPointShape();
1273    // resize the numeric array to the shape just calculated    size_t length=DataTypes::noValues(dataPointShape);
   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]);  
   }  
1274    
1275    // added for the MPI communication    // added for the MPI communication
   length=1;  
   for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];  
1276    double *tmpData = new double[length];    double *tmpData = new double[length];
1277    
1278    //    // updated for the MPI case
1279    // load the values for the data point into the numeric array.    if( get_MPIRank()==procNo ){
1280          if (getNumDataPointsPerSample()>0) {
1281      // updated for the MPI case      int sampleNo = dataPointNo/getNumDataPointsPerSample();
1282      if( get_MPIRank()==procNo ){      int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1283               if (getNumDataPointsPerSample()>0) {      //
1284                  int sampleNo = dataPointNo/getNumDataPointsPerSample();      // Check a valid sample number has been supplied
1285                  int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();      if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1286                  //          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;  
         }  
             }  
1287      }      }
1288          #ifdef PASO_MPI  
1289          // broadcast the data to all other processes      //
1290      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );      // Check a valid data point number has been supplied
1291          #endif      if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1292            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;  
1293      }      }
1294        // TODO: global error handling
1295        DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1296    
1297      delete [] tmpData;      memcpy(tmpData,&(getDataAtOffsetRO(offset)),length*sizeof(double));
1298         }
1299      }
1300    #ifdef PASO_MPI
1301      // broadcast the data to all other processes
1302      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1303    #endif
1304    
1305      boost::python::tuple t=pointToTuple(dataPointShape,tmpData);
1306      delete [] tmpData;
1307    //    //
1308    // return the loaded array    // return the loaded array
1309    return numArray;    return t;
1310    
1311  }  }
1312    
1313    
1314  boost::python::numeric::array  boost::python::object
1315  Data::integrate_const() const  Data::integrateToTuple_const() const
1316  {  {
1317    if (isLazy())    if (isLazy())
1318    {    {
# Line 1127  Data::integrate_const() const Line 1321  Data::integrate_const() const
1321    return integrateWorker();    return integrateWorker();
1322  }  }
1323    
1324  boost::python::numeric::array  boost::python::object
1325  Data::integrate()  Data::integrateToTuple()
1326  {  {
1327    if (isLazy())    if (isLazy())
1328    {    {
1329      expand();      expand();   // Can't do a non-resolving version of this without changing the domain object
1330    }    }         // see the dom->setToIntegrals call. Not saying it can't be done, just not doing it yet.
1331    return integrateWorker();    return integrateWorker();
 }  
   
1332    
1333    }
1334    
1335  boost::python::numeric::array  boost::python::object
1336  Data::integrateWorker() const  Data::integrateWorker() const
1337  {  {
   int index;  
   int rank = getDataPointRank();  
1338    DataTypes::ShapeType shape = getDataPointShape();    DataTypes::ShapeType shape = getDataPointShape();
1339    int dataPointSize = getDataPointSize();    int dataPointSize = getDataPointSize();
1340    
# Line 1151  Data::integrateWorker() const Line 1342  Data::integrateWorker() const
1342    // calculate the integral values    // calculate the integral values
1343    vector<double> integrals(dataPointSize);    vector<double> integrals(dataPointSize);
1344    vector<double> integrals_local(dataPointSize);    vector<double> integrals_local(dataPointSize);
1345      const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1346      if (dom==0)
1347      {            
1348        throw DataException("Can not integrate over non-continuous domains.");
1349      }
1350  #ifdef PASO_MPI  #ifdef PASO_MPI
1351    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals_local,*this);    dom->setToIntegrals(integrals_local,*this);
1352    // 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
1353    double *tmp = new double[dataPointSize];    double *tmp = new double[dataPointSize];
1354    double *tmp_local = new double[dataPointSize];    double *tmp_local = new double[dataPointSize];
1355    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]; }
1356    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 );
1357    for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }    for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1358      tuple result=pointToTuple(shape,tmp);
1359    delete[] tmp;    delete[] tmp;
1360    delete[] tmp_local;    delete[] tmp_local;
1361  #else  #else
1362    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);    dom->setToIntegrals(integrals,*this);
1363    /*  double *tmp = new double[dataPointSize];
1364      for (int i=0; i<dataPointSize; i++) { tmp[i]=integrals[i]; }*/
1365      tuple result=pointToTuple(shape,integrals);
1366    //   delete tmp;
1367  #endif  #endif
1368    
   //  
   // 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];  
           }  
         }  
       }  
     }  
   }  
1369    
1370    //    return result;
   // return the loaded array  
   return bp_array;  
1371  }  }
1372    
1373  Data  Data
1374  Data::sin() const  Data::sin() const
1375  {  {
1376    if (isLazy())    MAKELAZYOP(SIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIN);  
     return Data(c);  
   }  
1377    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1378  }  }
1379    
1380  Data  Data
1381  Data::cos() const  Data::cos() const
1382  {  {
1383    if (isLazy())    MAKELAZYOP(COS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COS);  
     return Data(c);  
   }  
1384    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1385  }  }
1386    
1387  Data  Data
1388  Data::tan() const  Data::tan() const
1389  {  {
1390    if (isLazy())    MAKELAZYOP(TAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TAN);  
     return Data(c);  
   }  
1391    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1392  }  }
1393    
1394  Data  Data
1395  Data::asin() const  Data::asin() const
1396  {  {
1397    if (isLazy())    MAKELAZYOP(ASIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ASIN);  
     return Data(c);  
   }  
1398    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1399  }  }
1400    
1401  Data  Data
1402  Data::acos() const  Data::acos() const
1403  {  {
1404    if (isLazy())    MAKELAZYOP(ACOS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOS);  
     return Data(c);  
   }  
1405    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1406  }  }
1407    
# Line 1279  Data::acos() const Line 1409  Data::acos() const
1409  Data  Data
1410  Data::atan() const  Data::atan() const
1411  {  {
1412    if (isLazy())    MAKELAZYOP(ATAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ATAN);  
     return Data(c);  
   }  
1413    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1414  }  }
1415    
1416  Data  Data
1417  Data::sinh() const  Data::sinh() const
1418  {  {
1419    if (isLazy())    MAKELAZYOP(SINH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SINH);  
     return Data(c);  
   }  
1420    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1421  }  }
1422    
1423  Data  Data
1424  Data::cosh() const  Data::cosh() const
1425  {  {
1426    if (isLazy())    MAKELAZYOP(COSH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COSH);  
     return Data(c);  
   }  
1427    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1428  }  }
1429    
1430  Data  Data
1431  Data::tanh() const  Data::tanh() const
1432  {  {
1433    if (isLazy())    MAKELAZYOP(TANH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TANH);  
     return Data(c);  
   }  
1434    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1435  }  }
1436    
# Line 1327  Data::erf() const Line 1441  Data::erf() const
1441  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1442    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1443  #else  #else
1444    if (isLazy())    MAKELAZYOP(ERF)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ERF);  
     return Data(c);  
   }  
1445    return C_TensorUnaryOperation(*this, ::erf);    return C_TensorUnaryOperation(*this, ::erf);
1446  #endif  #endif
1447  }  }
# Line 1339  Data::erf() const Line 1449  Data::erf() const
1449  Data  Data
1450  Data::asinh() const  Data::asinh() const
1451  {  {
1452    if (isLazy())    MAKELAZYOP(ASINH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);  
     return Data(c);  
   }  
1453  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1454    return C_TensorUnaryOperation(*this, escript::asinh_substitute);    return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1455  #else  #else
# Line 1354  Data::asinh() const Line 1460  Data::asinh() const
1460  Data  Data
1461  Data::acosh() const  Data::acosh() const
1462  {  {
1463    if (isLazy())    MAKELAZYOP(ACOSH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);  
     return Data(c);  
   }  
1464  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1465    return C_TensorUnaryOperation(*this, escript::acosh_substitute);    return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1466  #else  #else
# Line 1369  Data::acosh() const Line 1471  Data::acosh() const
1471  Data  Data
1472  Data::atanh() const  Data::atanh() const
1473  {  {
1474    if (isLazy())    MAKELAZYOP(ATANH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);  
     return Data(c);  
   }  
1475  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1476    return C_TensorUnaryOperation(*this, escript::atanh_substitute);    return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1477  #else  #else
# Line 1383  Data::atanh() const Line 1481  Data::atanh() const
1481    
1482  Data  Data
1483  Data::log10() const  Data::log10() const
1484  {  if (isLazy())  {
1485    {    MAKELAZYOP(LOG10)
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG10);  
     return Data(c);  
   }  
1486    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1487  }  }
1488    
1489  Data  Data
1490  Data::log() const  Data::log() const
1491  {  {
1492    if (isLazy())    MAKELAZYOP(LOG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG);  
     return Data(c);  
   }  
1493    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1494  }  }
1495    
1496  Data  Data
1497  Data::sign() const  Data::sign() const
1498  {  {
1499    if (isLazy())    MAKELAZYOP(SIGN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIGN);  
     return Data(c);  
   }  
1500    return C_TensorUnaryOperation(*this, escript::fsign);    return C_TensorUnaryOperation(*this, escript::fsign);
1501  }  }
1502    
1503  Data  Data
1504  Data::abs() const  Data::abs() const
1505  {  {
1506    if (isLazy())    MAKELAZYOP(ABS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ABS);  
     return Data(c);  
   }  
1507    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1508  }  }
1509    
1510  Data  Data
1511  Data::neg() const  Data::neg() const
1512  {  {
1513    if (isLazy())    MAKELAZYOP(NEG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NEG);  
     return Data(c);  
   }  
1514    return C_TensorUnaryOperation(*this, negate<double>());    return C_TensorUnaryOperation(*this, negate<double>());
1515  }  }
1516    
# Line 1448  Data::pos() const Line 1527  Data::pos() const
1527    
1528  Data  Data
1529  Data::exp() const  Data::exp() const
1530  {    {
1531    if (isLazy())    MAKELAZYOP(EXP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),EXP);  
     return Data(c);  
   }  
1532    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1533  }  }
1534    
1535  Data  Data
1536  Data::sqrt() const  Data::sqrt() const
1537  {  {
1538    if (isLazy())    MAKELAZYOP(SQRT)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SQRT);  
     return Data(c);  
   }  
1539    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1540  }  }
1541    
# Line 1483  Data::Lsup() Line 1554  Data::Lsup()
1554  {  {
1555     if (isLazy())     if (isLazy())
1556     {     {
1557      expand();      if (!actsExpanded() || CHECK_DO_CRES)
1558        {
1559            resolve();
1560        }
1561        else
1562        {
1563    #ifdef PASO_MPI
1564            return lazyAlgWorker<AbsMax>(0,MPI_MAX);
1565    #else
1566            return lazyAlgWorker<AbsMax>(0);
1567    #endif
1568        }
1569     }     }
1570     return LsupWorker();     return LsupWorker();
1571  }  }
# Line 1503  Data::sup() Line 1585  Data::sup()
1585  {  {
1586     if (isLazy())     if (isLazy())
1587     {     {
1588      expand();      if (!actsExpanded() || CHECK_DO_CRES)
1589        {
1590            resolve();
1591        }
1592        else
1593        {
1594    #ifdef PASO_MPI
1595            return lazyAlgWorker<FMax>(numeric_limits<double>::max()*-1, MPI_MAX);
1596    #else
1597            return lazyAlgWorker<FMax>(numeric_limits<double>::max()*-1);
1598    #endif
1599        }
1600     }     }
1601     return supWorker();     return supWorker();
1602  }  }
# Line 1523  Data::inf() Line 1616  Data::inf()
1616  {  {
1617     if (isLazy())     if (isLazy())
1618     {     {
1619      expand();      if (!actsExpanded() || CHECK_DO_CRES)
1620        {
1621            resolve();
1622        }
1623        else
1624        {
1625    #ifdef PASO_MPI
1626            return lazyAlgWorker<FMin>(numeric_limits<double>::max(), MPI_MIN);
1627    #else
1628            return lazyAlgWorker<FMin>(numeric_limits<double>::max());
1629    #endif
1630        }
1631     }     }
1632     return infWorker();     return infWorker();
1633  }  }
1634    
1635    template <class BinaryOp>
1636    double
1637    #ifdef PASO_MPI
1638    Data::lazyAlgWorker(double init, MPI_Op mpiop_type)
1639    #else
1640    Data::lazyAlgWorker(double init)
1641    #endif
1642    {
1643       if (!isLazy() || !m_data->actsExpanded())
1644       {
1645        throw DataException("Error - lazyAlgWorker can only be called on lazy(expanded) data.");
1646       }
1647       DataLazy* dl=dynamic_cast<DataLazy*>(m_data.get());
1648       EsysAssert((dl!=0), "Programming error - casting to DataLazy.");
1649       BufferGroup* bg=allocSampleBuffer();
1650       double val=init;
1651       int i=0;
1652       const size_t numsamples=getNumSamples();
1653       const size_t samplesize=getNoValues()*getNumDataPointsPerSample();
1654       BinaryOp operation;
1655       #pragma omp parallel private(i)
1656       {
1657        double localtot=init;
1658        #pragma omp for schedule(static) private(i)
1659        for (i=0;i<numsamples;++i)
1660        {
1661            size_t roffset=0;
1662            const DataTypes::ValueType* v=dl->resolveSample(*bg, i, roffset);
1663            // Now we have the sample, run operation on all points
1664            for (size_t j=0;j<samplesize;++j)
1665            {
1666            localtot=operation(localtot,(*v)[j+roffset]);
1667            }
1668        }
1669        #pragma omp critical
1670        val=operation(val,localtot);
1671       }
1672       freeSampleBuffer(bg);
1673    #ifdef PASO_MPI
1674       double globalValue;
1675       MPI_Allreduce( &val, &globalValue, 1, MPI_DOUBLE, mpiop_type, MPI_COMM_WORLD );
1676       return globalValue;
1677    #else
1678       return val;
1679    #endif
1680    }
1681    
1682  double  double
1683  Data::LsupWorker() const  Data::LsupWorker() const
1684  {  {
# Line 1580  Data::infWorker() const Line 1731  Data::infWorker() const
1731  #endif  #endif
1732  }  }
1733    
 /* TODO */  
1734  /* global reduction */  /* global reduction */
1735  Data  
1736  Data::maxval() const  
1737    inline Data
1738    Data::minval_nonlazy() const
1739    {
1740      //
1741      // set the initial minimum value to max possible double
1742      FMin fmin_func;
1743      return dp_algorithm(fmin_func,numeric_limits<double>::max());
1744    }
1745    
1746    
1747    inline Data
1748    Data::maxval_nonlazy() const
1749  {  {
   if (isLazy())  
   {  
     Data temp(*this);   // to get around the fact that you can't resolve a const Data  
     temp.resolve();  
     return temp.maxval();  
   }  
1750    //    //
1751    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1752    FMax fmax_func;    FMax fmax_func;
1753    return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);    return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1754  }  }
1755    
1756    
1757    
1758    Data
1759    Data::maxval() const
1760    {
1761       MAKELAZYOP(MAXVAL)
1762       return maxval_nonlazy();
1763    }
1764    
1765    
1766  Data  Data
1767  Data::minval() const  Data::minval() const
1768  {  {
1769    if (isLazy())    MAKELAZYOP(MINVAL)
1770    {    return minval_nonlazy();
     Data temp(*this);   // to get around the fact that you can't resolve a const Data  
     temp.resolve();  
     return temp.minval();  
   }  
   //  
   // set the initial minimum value to max possible double  
   FMin fmin_func;  
   return dp_algorithm(fmin_func,numeric_limits<double>::max());  
1771  }  }
1772    
1773    
1774  Data  Data
1775  Data::swapaxes(const int axis0, const int axis1) const  Data::swapaxes(const int axis0, const int axis1) const
1776  {  {
# Line 1633  Data::swapaxes(const int axis0, const in Line 1792  Data::swapaxes(const int axis0, const in
1792       if (axis0 == axis1) {       if (axis0 == axis1) {
1793           throw DataException("Error - Data::swapaxes: axis indices must be different.");           throw DataException("Error - Data::swapaxes: axis indices must be different.");
1794       }       }
1795       if (axis0 > axis1) {       MAKELAZYOP2(SWAP,axis0,axis1)
1796           axis0_tmp=axis1;       if (axis0 > axis1)
1797           axis1_tmp=axis0;       {
1798       } else {      axis0_tmp=axis1;
1799           axis0_tmp=axis0;      axis1_tmp=axis0;
          axis1_tmp=axis1;  
1800       }       }
1801       for (int i=0; i<rank; i++) {       else
1802         if (i == axis0_tmp) {       {
1803            ev_shape.push_back(s[axis1_tmp]);      axis0_tmp=axis0;
1804         } else if (i == axis1_tmp) {      axis1_tmp=axis1;
1805            ev_shape.push_back(s[axis0_tmp]);       }
1806         } else {       for (int i=0; i<rank; i++)
1807            ev_shape.push_back(s[i]);       {
1808         }      if (i == axis0_tmp)
1809        {
1810            ev_shape.push_back(s[axis1_tmp]);
1811        }
1812        else if (i == axis1_tmp)
1813        {
1814            ev_shape.push_back(s[axis0_tmp]);
1815        }
1816        else
1817        {
1818            ev_shape.push_back(s[i]);
1819        }
1820       }       }
1821       Data ev(0.,ev_shape,getFunctionSpace());       Data ev(0.,ev_shape,getFunctionSpace());
1822       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1823       m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);       m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1824       return ev;       return ev;
   
1825  }  }
1826    
1827  Data  Data
# Line 1672  Data::symmetric() const Line 1840  Data::symmetric() const
1840       else {       else {
1841          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.");
1842       }       }
1843       if (isLazy())       MAKELAZYOP(SYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SYM);  
     return Data(c);  
      }  
1844       Data ev(0.,getDataPointShape(),getFunctionSpace());       Data ev(0.,getDataPointShape(),getFunctionSpace());
1845       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1846       m_data->symmetric(ev.m_data.get());       m_data->symmetric(ev.m_data.get());
# Line 1686  Data::symmetric() const Line 1850  Data::symmetric() const
1850  Data  Data
1851  Data::nonsymmetric() const  Data::nonsymmetric() const
1852  {  {
1853       if (isLazy())       MAKELAZYOP(NSYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NSYM);  
     return Data(c);  
      }  
1854       // check input       // check input
1855       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1856       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1722  Data::nonsymmetric() const Line 1882  Data::nonsymmetric() const
1882       }       }
1883  }  }
1884    
   
 // 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).  
1885  Data  Data
1886  Data::trace(int axis_offset) const  Data::trace(int axis_offset) const
1887  {  {    
1888       if (isLazy())       MAKELAZYOPOFF(TRACE,axis_offset)
1889         if ((axis_offset<0) || (axis_offset>getDataPointRank()))
1890       {       {
1891      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);  
1892       }       }
1893       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1894       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1787  Data::trace(int axis_offset) const Line 1940  Data::trace(int axis_offset) const
1940  Data  Data
1941  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1942  {      {    
1943       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);  
      }  
1944       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1945       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
1946       // 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 1802  Data::transpose(int axis_offset) const Line 1950  Data::transpose(int axis_offset) const
1950          throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);          throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1951       }       }
1952       for (int i=0; i<rank; i++) {       for (int i=0; i<rank; i++) {
1953    
1954         int index = (axis_offset+i)%rank;         int index = (axis_offset+i)%rank;
1955         ev_shape.push_back(s[index]); // Append to new shape         ev_shape.push_back(s[index]); // Append to new shape
1956       }       }
# Line 1886  Data::calc_minGlobalDataPoint(int& ProcN Line 2035  Data::calc_minGlobalDataPoint(int& ProcN
2035    int lowi=0,lowj=0;    int lowi=0,lowj=0;
2036    double min=numeric_limits<double>::max();    double min=numeric_limits<double>::max();
2037    
2038    Data temp=minval();    Data temp=minval_nonlazy();   // need to do this to prevent autolazy from reintroducing laziness
2039    
2040    int numSamples=temp.getNumSamples();    int numSamples=temp.getNumSamples();
2041    int numDPPSample=temp.getNumDataPointsPerSample();    int numDPPSample=temp.getNumDataPointsPerSample();
# Line 1894  Data::calc_minGlobalDataPoint(int& ProcN Line 2043  Data::calc_minGlobalDataPoint(int& ProcN
2043    double next,local_min;    double next,local_min;
2044    int local_lowi=0,local_lowj=0;        int local_lowi=0,local_lowj=0;    
2045    
2046    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)    #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
2047    {    {
2048      local_min=min;      local_min=min;
2049      #pragma omp for private(i,j) schedule(static)      #pragma omp for private(i,j) schedule(static)
2050      for (i=0; i<numSamples; i++) {      for (i=0; i<numSamples; i++) {
2051        for (j=0; j<numDPPSample; j++) {        for (j=0; j<numDPPSample; j++) {
2052          next=temp.getDataAtOffset(temp.getDataOffset(i,j));          next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2053          if (next<local_min) {          if (next<local_min) {
2054            local_min=next;            local_min=next;
2055            local_lowi=i;            local_lowi=i;
# Line 1909  Data::calc_minGlobalDataPoint(int& ProcN Line 2058  Data::calc_minGlobalDataPoint(int& ProcN
2058        }        }
2059      }      }
2060      #pragma omp critical      #pragma omp critical
2061      if (local_min<min) {      if (local_min<min) {    // If we found a smaller value than our sentinel
2062        min=local_min;        min=local_min;
2063        lowi=local_lowi;        lowi=local_lowi;
2064        lowj=local_lowj;        lowj=local_lowj;
# Line 1917  Data::calc_minGlobalDataPoint(int& ProcN Line 2066  Data::calc_minGlobalDataPoint(int& ProcN
2066    }    }
2067    
2068  #ifdef PASO_MPI  #ifdef PASO_MPI
2069      // determine the processor on which the minimum occurs    // determine the processor on which the minimum occurs
2070      next = temp.getDataPoint(lowi,lowj);    next = temp.getDataPointRO(lowi,lowj);
2071      int lowProc = 0;    int lowProc = 0;
2072      double *globalMins = new double[get_MPISize()+1];    double *globalMins = new double[get_MPISize()+1];
2073      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );    int error;
2074      error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
2075      if( get_MPIRank()==0 ){  
2076          next = globalMins[lowProc];    if( get_MPIRank()==0 ){
2077          for( i=1; i<get_MPISize(); i++ )      next = globalMins[lowProc];
2078              if( next>globalMins[i] ){      for( i=1; i<get_MPISize(); i++ )
2079                  lowProc = i;          if( next>globalMins[i] ){
2080                  next = globalMins[i];              lowProc = i;
2081              }              next = globalMins[i];
2082            }
2083      }
2084      MPI_Bcast( &lowProc, 1, MPI_INT, 0, get_MPIComm() );
2085      DataPointNo = lowj + lowi * numDPPSample;
2086      MPI_Bcast(&DataPointNo, 1, MPI_INT, lowProc, get_MPIComm() );
2087      delete [] globalMins;
2088      ProcNo = lowProc;
2089    #else
2090      ProcNo = 0;
2091      DataPointNo = lowj + lowi * numDPPSample;
2092    #endif
2093    }
2094    
2095    
2096    const boost::python::tuple
2097    Data::maxGlobalDataPoint() const
2098    {
2099      int DataPointNo;
2100      int ProcNo;
2101      calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2102      return make_tuple(ProcNo,DataPointNo);
2103    }
2104    
2105    void
2106    Data::calc_maxGlobalDataPoint(int& ProcNo,
2107                            int& DataPointNo) const
2108    {
2109      if (isLazy())
2110      {
2111        Data temp(*this);   // to get around the fact that you can't resolve a const Data
2112        temp.resolve();
2113        return temp.calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2114      }
2115      int i,j;
2116      int highi=0,highj=0;
2117    //-------------
2118      double max= -numeric_limits<double>::max();
2119    
2120      Data temp=maxval_nonlazy();   // need to do this to prevent autolazy from reintroducing laziness
2121    
2122      int numSamples=temp.getNumSamples();
2123      int numDPPSample=temp.getNumDataPointsPerSample();
2124    
2125      double next,local_max;
2126      int local_highi=0,local_highj=0;  
2127    
2128      #pragma omp parallel firstprivate(local_highi,local_highj) private(next,local_max)
2129      {
2130        local_max=max;
2131        #pragma omp for private(i,j) schedule(static)
2132        for (i=0; i<numSamples; i++) {
2133          for (j=0; j<numDPPSample; j++) {
2134            next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2135            if (next>local_max) {
2136              local_max=next;
2137              local_highi=i;
2138              local_highj=j;
2139            }
2140          }
2141        }
2142        #pragma omp critical
2143        if (local_max>max) {    // If we found a larger value than our sentinel
2144          max=local_max;
2145          highi=local_highi;
2146          highj=local_highj;
2147        }
2148      }
2149    #ifdef PASO_MPI
2150      // determine the processor on which the maximum occurs
2151      next = temp.getDataPointRO(highi,highj);
2152      int highProc = 0;
2153      double *globalMaxs = new double[get_MPISize()+1];
2154      int error;
2155      error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMaxs, 1, MPI_DOUBLE, 0, get_MPIComm() );
2156      if( get_MPIRank()==0 ){
2157        next = globalMaxs[highProc];
2158        for( i=1; i<get_MPISize(); i++ )
2159        {
2160        if( next<globalMaxs[i] )
2161        {
2162            highProc = i;
2163            next = globalMaxs[i];
2164      }      }
2165      MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );      }
2166      }
2167      MPI_Bcast( &highProc, 1, MPI_INT, 0, get_MPIComm() );
2168      DataPointNo = highj + highi * numDPPSample;  
2169      MPI_Bcast(&DataPointNo, 1, MPI_INT, highProc, get_MPIComm() );
2170    
2171      delete [] globalMins;    delete [] globalMaxs;
2172      ProcNo = lowProc;    ProcNo = highProc;
2173  #else  #else
2174      ProcNo = 0;    ProcNo = 0;
2175      DataPointNo = highj + highi * numDPPSample;
2176  #endif  #endif
   DataPointNo = lowj + lowi * numDPPSample;  
2177  }  }
2178    
2179  void  void
# Line 1977  Data::saveVTK(std::string fileName) cons Line 2212  Data::saveVTK(std::string fileName) cons
2212    }    }
2213    boost::python::dict args;    boost::python::dict args;
2214    args["data"]=boost::python::object(this);    args["data"]=boost::python::object(this);
2215    getDomain()->saveVTK(fileName,args);    getDomain()->saveVTK(fileName,args,"","");
2216    return;    return;
2217  }  }
2218    
2219    
2220    
2221  Data&  Data&
2222  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
2223  {  {
2224    if (isProtected()) {    if (isProtected()) {
2225          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2226    }    }
2227    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,ADD)    // for lazy + is equivalent to +=
2228    {    exclusiveWrite();         // Since Lazy data does not modify its leaves we only need to worry here
2229      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),ADD); // for lazy + is equivalent to +=    binaryOp(right,plus<double>());
2230          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,plus<double>());  
     return (*this);  
   }  
2231  }  }
2232    
2233  Data&  Data&
# Line 2007  Data::operator+=(const boost::python::ob Line 2237  Data::operator+=(const boost::python::ob
2237          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2238    }    }
2239    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2240    if (isLazy())    (*this)+=tmp;
2241    {    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);  
   }  
2242  }  }
2243    
2244  // Hmmm, operator= makes a deep copy but the copy constructor does not?  // Hmmm, operator= makes a deep copy but the copy constructor does not?
2245  Data&  Data&
2246  Data::operator=(const Data& other)  Data::operator=(const Data& other)
2247  {  {
2248    copy(other);    m_protected=false;        // since any changes should be caught by exclusiveWrite();
2249    //   m_data=other.m_data;
2250      set_m_data(other.m_data);
2251    return (*this);    return (*this);
2252  }  }
2253    
# Line 2034  Data::operator-=(const Data& right) Line 2257  Data::operator-=(const Data& right)
2257    if (isProtected()) {    if (isProtected()) {
2258          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2259    }    }
2260    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,SUB)
2261    {    exclusiveWrite();
2262      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),SUB); // for lazy - is equivalent to -=    binaryOp(right,minus<double>());
2263          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,minus<double>());  
     return (*this);  
   }  
2264  }  }
2265    
2266  Data&  Data&
# Line 2054  Data::operator-=(const boost::python::ob Line 2270  Data::operator-=(const boost::python::ob
2270          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2271    }    }
2272    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2273    if (isLazy())    (*this)-=tmp;
2274    {    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);  
   }  
2275  }  }
2276    
2277  Data&  Data&
# Line 2073  Data::operator*=(const Data& right) Line 2280  Data::operator*=(const Data& right)
2280    if (isProtected()) {    if (isProtected()) {
2281          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2282    }    }
2283    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,MUL)
2284    {    exclusiveWrite();
2285      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),MUL); // for lazy * is equivalent to *=    binaryOp(right,multiplies<double>());
2286          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,multiplies<double>());  
     return (*this);  
   }  
2287  }  }
2288    
2289  Data&  Data&
# Line 2093  Data::operator*=(const boost::python::ob Line 2293  Data::operator*=(const boost::python::ob
2293          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2294    }    }
2295    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2296    if (isLazy())    (*this)*=tmp;
2297    {    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);  
   }  
2298  }  }
2299    
2300  Data&  Data&
# Line 2112  Data::operator/=(const Data& right) Line 2303  Data::operator/=(const Data& right)
2303    if (isProtected()) {    if (isProtected()) {
2304          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2305    }    }
2306    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,DIV)
2307    {    exclusiveWrite();
2308      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),DIV); // for lazy / is equivalent to /=    binaryOp(right,divides<double>());
2309          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,divides<double>());  
     return (*this);  
   }  
2310  }  }
2311    
2312  Data&  Data&
# Line 2132  Data::operator/=(const boost::python::ob Line 2316  Data::operator/=(const boost::python::ob
2316          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2317    }    }
2318    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2319    if (isLazy())    (*this)/=tmp;
2320    {    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);  
   }  
2321  }  }
2322    
2323  Data  Data
# Line 2162  Data::powO(const boost::python::object& Line 2337  Data::powO(const boost::python::object&
2337  Data  Data
2338  Data::powD(const Data& right) const  Data::powD(const Data& right) const
2339  {  {
2340    if (isLazy() || right.isLazy())    MAKELAZYBIN(right,POW)
   {  
     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),POW);  
     return Data(c);  
   }  
2341    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2342  }  }
2343    
# Line 2175  Data::powD(const Data& right) const Line 2346  Data::powD(const Data& right) const
2346  Data  Data
2347  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
2348  {  {
2349    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,ADD)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);  
     return Data(c);  
   }  
2350    return C_TensorBinaryOperation(left, right, plus<double>());    return C_TensorBinaryOperation(left, right, plus<double>());
2351  }  }
2352    
# Line 2188  escript::operator+(const Data& left, con Line 2355  escript::operator+(const Data& left, con
2355  Data  Data
2356  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
2357  {  {
2358    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,SUB)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);  
     return Data(c);  
   }  
2359    return C_TensorBinaryOperation(left, right, minus<double>());    return C_TensorBinaryOperation(left, right, minus<double>());
2360  }  }
2361    
# Line 2201  escript::operator-(const Data& left, con Line 2364  escript::operator-(const Data& left, con
2364  Data  Data
2365  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
2366  {  {
2367    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,MUL)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);  
     return Data(c);  
   }  
2368    return C_TensorBinaryOperation(left, right, multiplies<double>());    return C_TensorBinaryOperation(left, right, multiplies<double>());
2369  }  }
2370    
# Line 2214  escript::operator*(const Data& left, con Line 2373  escript::operator*(const Data& left, con
2373  Data  Data
2374  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
2375  {  {
2376    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,DIV)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),DIV);  
     return Data(c);  
   }  
2377    return C_TensorBinaryOperation(left, right, divides<double>());    return C_TensorBinaryOperation(left, right, divides<double>());
2378  }  }
2379    
# Line 2227  escript::operator/(const Data& left, con Line 2382  escript::operator/(const Data& left, con
2382  Data  Data
2383  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
2384  {  {
2385    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2386    {    MAKELAZYBIN2(left,tmp,ADD)
2387      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);  
2388  }  }
2389    
2390  //  //
# Line 2240  escript::operator+(const Data& left, con Line 2392  escript::operator+(const Data& left, con
2392  Data  Data
2393  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
2394  {  {
2395    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2396    {    MAKELAZYBIN2(left,tmp,SUB)
2397      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);  
2398  }  }
2399    
2400  //  //
# Line 2253  escript::operator-(const Data& left, con Line 2402  escript::operator-(const Data& left, con
2402  Data  Data
2403  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
2404  {  {
2405    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2406    {    MAKELAZYBIN2(left,tmp,MUL)
2407      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);  
2408  }  }
2409    
2410  //  //
# Line 2266  escript::operator*(const Data& left, con Line 2412  escript::operator*(const Data& left, con
2412  Data  Data
2413  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
2414  {  {
2415    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2416    {    MAKELAZYBIN2(left,tmp,DIV)
2417      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);  
2418  }  }
2419    
2420  //  //
# Line 2279  escript::operator/(const Data& left, con Line 2422  escript::operator/(const Data& left, con
2422  Data  Data
2423  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
2424  {  {
2425    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2426    {    MAKELAZYBIN2(tmp,right,ADD)
2427      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;  
2428  }  }
2429    
2430  //  //
# Line 2292  escript::operator+(const boost::python:: Line 2432  escript::operator+(const boost::python::
2432  Data  Data
2433  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
2434  {  {
2435    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2436    {    MAKELAZYBIN2(tmp,right,SUB)
2437      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;  
2438  }  }
2439    
2440  //  //
# Line 2305  escript::operator-(const boost::python:: Line 2442  escript::operator-(const boost::python::
2442  Data  Data
2443  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
2444  {  {
2445    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2446    {    MAKELAZYBIN2(tmp,right,MUL)
2447      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;  
2448  }  }
2449    
2450  //  //
# Line 2318  escript::operator*(const boost::python:: Line 2452  escript::operator*(const boost::python::
2452  Data  Data
2453  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
2454  {  {
2455    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2456    {    MAKELAZYBIN2(tmp,right,DIV)
2457      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;  
2458  }  }
2459    
2460    
# Line 2364  void Line 2495  void
2495  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
2496                 const Data& value)                 const Data& value)
2497  {  {
 //  const DataArrayView& view=getPointDataView();  
   
2498    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2499    if (slice_region.size()!=getDataPointRank()) {    if (slice_region.size()!=getDataPointRank()) {
2500      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
2501    }    }
2502      exclusiveWrite();
2503    if (getFunctionSpace()!=value.getFunctionSpace()) {    if (getFunctionSpace()!=value.getFunctionSpace()) {
2504       setSlice(Data(value,getFunctionSpace()),slice_region);       setSlice(Data(value,getFunctionSpace()),slice_region);
2505    } else {    } else {
# Line 2384  Data::setSlice(const Data& value, Line 2514  Data::setSlice(const Data& value,
2514    if (isProtected()) {    if (isProtected()) {
2515          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2516    }    }
2517    FORCERESOLVE;    forceResolve();
2518  /*  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.");  
   }*/  
2519    Data tempValue(value);    Data tempValue(value);
2520    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
2521    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 2431  Data::typeMatchRight(const Data& right) Line 2558  Data::typeMatchRight(const Data& right)
2558    }    }
2559  }  }
2560    
2561    // The normal TaggedValue adds the tag if it is not already present
2562    // This form does not. It throws instead.
2563    // This is because the names are maintained by the domain and cannot be added
2564    // without knowing the tag number to map it to.
2565  void  void
2566  Data::setTaggedValueByName(std::string name,  Data::setTaggedValueByName(std::string name,
2567                             const boost::python::object& value)                             const boost::python::object& value)
2568  {  {
2569       if (getFunctionSpace().getDomain()->isValidTagName(name)) {       if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2570      FORCERESOLVE;      forceResolve();
2571        exclusiveWrite();
2572          int tagKey=getFunctionSpace().getDomain()->getTag(name);          int tagKey=getFunctionSpace().getDomain()->getTag(name);
2573          setTaggedValue(tagKey,value);          setTaggedValue(tagKey,value);
2574       }       }
2575         else
2576         {                  // The
2577        throw DataException("Error - unknown tag in setTaggedValueByName.");
2578         }
2579  }  }
2580    
2581  void  void
2582  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
2583                       const boost::python::object& value)                       const boost::python::object& value)
# Line 2450  Data::setTaggedValue(int tagKey, Line 2587  Data::setTaggedValue(int tagKey,
2587    }    }
2588    //    //
2589    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2590    FORCERESOLVE;    forceResolve();
2591      exclusiveWrite();
2592    if (isConstant()) tag();    if (isConstant()) tag();
2593    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]));  
   }  
2594    
2595    DataVector temp_data2;    DataVector temp_data2;
2596    temp_data2.copyFromNumArray(asNumArray);    temp_data2.copyFromArray(w,1);
2597    
2598    m_data->setTaggedValue(tagKey,tempShape, temp_data2);    m_data->setTaggedValue(tagKey,w.getShape(), temp_data2);
2599  }  }
2600    
2601    
# Line 2478  Data::setTaggedValueFromCPP(int tagKey, Line 2610  Data::setTaggedValueFromCPP(int tagKey,
2610    }    }
2611    //    //
2612    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2613    FORCERESOLVE;    forceResolve();
2614    if (isConstant()) tag();    if (isConstant()) tag();
2615      exclusiveWrite();
2616    //    //
2617    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
2618    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
# Line 2512  escript::C_GeneralTensorProduct(Data& ar Line 2645  escript::C_GeneralTensorProduct(Data& ar
2645    // 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().
2646    
2647    // deal with any lazy data    // deal with any lazy data
2648    if (arg_0.isLazy()) {arg_0.resolve();}  //   if (arg_0.isLazy()) {arg_0.resolve();}
2649    if (arg_1.isLazy()) {arg_1.resolve();}  //   if (arg_1.isLazy()) {arg_1.resolve();}
2650      if (arg_0.isLazy() || arg_1.isLazy() || (AUTOLAZYON && (arg_0.isExpanded() || arg_1.isExpanded())))
2651      {
2652        DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2653        return Data(c);
2654      }
2655    
2656    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2657    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
# Line 2590  escript::C_GeneralTensorProduct(Data& ar Line 2728  escript::C_GeneralTensorProduct(Data& ar
2728       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
2729    }    }
2730    
2731      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2732      {
2733         ostringstream os;
2734         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2735         throw DataException(os.str());
2736      }
2737    
2738    // Declare output Data object    // Declare output Data object
2739    Data res;    Data res;
2740    
2741    if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {    if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2742      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2743      double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2744      double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2745      double *ptr_2 = &(res.getDataAtOffset(0));      double *ptr_2 = &(res.getDataAtOffsetRW(0));
2746      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);
2747    }    }
2748    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 2763  escript::C_GeneralTensorProduct(Data& ar
2763    
2764      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2765      int offset_0 = tmp_0->getPointOffset(0,0);      int offset_0 = tmp_0->getPointOffset(0,0);
2766      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));  
2767    
2768        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2769        double *ptr_2 = &(tmp_2->getDefaultValueRW(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 2637  escript::C_GeneralTensorProduct(Data& ar Line 2775  escript::C_GeneralTensorProduct(Data& ar
2775      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
2776      for (i=lookup_1.begin();i!=lookup_1.end();i++) {      for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2777        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);  
2778    
2779        double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2780        double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2781            
2782        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);
2783      }      }
# Line 2667  escript::C_GeneralTensorProduct(Data& ar Line 2801  escript::C_GeneralTensorProduct(Data& ar
2801        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2802          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2803          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2804          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2805          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2806          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2807          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);
2808        }        }
2809      }      }
# Line 2693  escript::C_GeneralTensorProduct(Data& ar Line 2827  escript::C_GeneralTensorProduct(Data& ar
2827    
2828      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2829      int offset_1 = tmp_1->getPointOffset(0,0);      int offset_1 = tmp_1->getPointOffset(0,0);
2830      double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2831      // Get the views      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2832  //     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));  
2833    
2834      // Compute an MVP for the default      // Compute an MVP for the default
2835      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 2837  escript::C_GeneralTensorProduct(Data& ar
2837      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2838      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
2839      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);  
2840    
2841        tmp_2->addTag(i->first);        tmp_2->addTag(i->first);
2842        double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2843        double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2844        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);
2845      }      }
2846    
# Line 2739  escript::C_GeneralTensorProduct(Data& ar Line 2861  escript::C_GeneralTensorProduct(Data& ar
2861      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2862      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2863    
2864  //     // Get the views      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2865  //     DataArrayView view_0 = tmp_0->getDefaultValue();      const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2866  //     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));  
   
2867    
2868      // Compute an MVP for the default      // Compute an MVP for the default
2869      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 2880  escript::C_GeneralTensorProduct(Data& ar
2880      // Compute an MVP for each tag      // Compute an MVP for each tag
2881      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2882      for (i=lookup_2.begin();i!=lookup_2.end();i++) {      for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2883  //       DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2884  //       DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2885  //       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));  
2886    
2887        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);
2888      }      }
# Line 2799  escript::C_GeneralTensorProduct(Data& ar Line 2904  escript::C_GeneralTensorProduct(Data& ar
2904      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2905      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2906        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
2907        double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2908        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2909          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2910          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2911          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2912          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2913          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);
2914        }        }
2915      }      }
# Line 2828  escript::C_GeneralTensorProduct(Data& ar Line 2933  escript::C_GeneralTensorProduct(Data& ar
2933        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2934          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2935          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2936          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2937          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2938          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2939          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);
2940        }        }
2941      }      }
# Line 2853  escript::C_GeneralTensorProduct(Data& ar Line 2958  escript::C_GeneralTensorProduct(Data& ar
2958      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2959      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2960        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2961        double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2962        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2963          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2964          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2965          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2966          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2967          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);
2968        }        }
2969      }      }
# Line 2883  escript::C_GeneralTensorProduct(Data& ar Line 2988  escript::C_GeneralTensorProduct(Data& ar
2988          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2989          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2990          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2991          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2992          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2993          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2994          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);
2995        }        }
2996      }      }
# Line 2923  Data::borrowReadyPtr() const Line 3028  Data::borrowReadyPtr() const
3028  std::string  std::string
3029  Data::toString() const  Data::toString() const
3030  {  {
3031      if (!m_data->isEmpty() &&      if (!m_data->isEmpty() &&
3032      getNumDataPoints()*getDataPointSize()>escriptParams.getInt("TOO_MANY_LINES"))      !m_data->isLazy() &&
3033        getLength()>escriptParams.getInt("TOO_MANY_LINES"))
3034      {      {
3035      stringstream temp;      stringstream temp;
3036      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 3040  Data::toString() const
3040  }  }
3041    
3042    
3043    // This method is not thread-safe
3044    DataTypes::ValueType::reference
3045    Data::getDataAtOffsetRW(DataTypes::ValueType::size_type i)
3046    {
3047        checkExclusiveWrite();
3048        return getReady()->getDataAtOffsetRW(i);
3049    }
3050    
3051    // This method is not thread-safe
3052  DataTypes::ValueType::const_reference  DataTypes::ValueType::const_reference
3053  Data::getDataAtOffset(DataTypes::ValueType::size_type i) const  Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i)
3054  {  {
3055      if (isLazy())      forceResolve();
3056      {      return getReady()->getDataAtOffsetRO(i);
     throw DataException("Programmer error - getDataAtOffset not permitted on lazy data (object is const which prevents resolving).");  
     }  
     return getReady()->getDataAtOffset(i);  
3057  }  }
3058    
3059    
3060  DataTypes::ValueType::reference  // DataTypes::ValueType::const_reference
3061  Data::getDataAtOffset(DataTypes::ValueType::size_type i)  // Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i) const
3062  {  // {
3063  //     if (isLazy())  //     if (isLazy())
3064  //     {  //     {
3065  //  throw DataException("getDataAtOffset not permitted on lazy data.");  //  throw DataException("Programmer error - getDataAtOffsetRO() not permitted on Lazy Data (object is const which prevents resolving).");
3066  //     }  //     }
3067      FORCERESOLVE;  //     return getReady()->getDataAtOffsetRO(i);
3068      return getReady()->getDataAtOffset(i);  // }
3069  }  
3070    
3071  DataTypes::ValueType::const_reference  DataTypes::ValueType::const_reference
3072  Data::getDataPoint(int sampleNo, int dataPointNo) const  Data::getDataPointRO(int sampleNo, int dataPointNo)
3073  {  {
3074      forceResolve();
3075    if (!isReady())    if (!isReady())
3076    {    {
3077      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.");
3078    }    }
3079    else    else
3080    {    {
3081      const DataReady* dr=getReady();      const DataReady* dr=getReady();
3082      return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));      return dr->getDataAtOffsetRO(dr->getPointOffset(sampleNo, dataPointNo));
3083    }    }
3084  }  }
3085    
3086    
3087    
3088    
3089  DataTypes::ValueType::reference  DataTypes::ValueType::reference
3090  Data::getDataPoint(int sampleNo, int dataPointNo)  Data::getDataPointRW(int sampleNo, int dataPointNo)
3091  {  {
3092    FORCERESOLVE;    checkExclusiveWrite();
3093    if (!isReady())    DataReady* dr=getReady();
3094    {    return dr->getDataAtOffsetRW(dr->getPointOffset(sampleNo, dataPointNo));
3095      throw DataException("Programmer error - getDataPoint() not permitted on Lazy Data.");  }
3096    }  
3097    else  BufferGroup*
3098    {  Data::allocSampleBuffer() const
3099      DataReady* dr=getReady();  {
3100      return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));       if (isLazy())
3101    }       {
3102        #ifdef _OPENMP
3103        int tnum=omp_get_max_threads();
3104        #else
3105        int tnum=1;
3106        #endif
3107        return new BufferGroup(getSampleBufferSize(),tnum);
3108         }
3109         else
3110         {
3111        return NULL;
3112         }
3113    }
3114    
3115    void
3116    Data::freeSampleBuffer(BufferGroup* bufferg)
3117    {
3118         if (bufferg!=0)
3119         {
3120        delete bufferg;
3121         }
3122    }
3123    
3124    
3125    Data
3126    Data::interpolateFromTable2DP(boost::python::object table, double Amin, double Astep,
3127            Data& B, double Bmin, double Bstep, double undef, bool check_boundaries)
3128    {
3129        WrappedArray t(table);
3130        return interpolateFromTable2D(t, Amin, Astep, undef, B, Bmin, Bstep,check_boundaries);
3131    }
3132    
3133    Data
3134    Data::interpolateFromTable1DP(boost::python::object table, double Amin, double Astep,
3135                              double undef,bool check_boundaries)
3136    {
3137        WrappedArray t(table);
3138        return interpolateFromTable1D(t, Amin, Astep, undef, check_boundaries);
3139    }
3140    
3141    
3142    Data
3143    Data::interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,
3144                     double undef, bool check_boundaries)
3145    {
3146        table.convertArray();       // critical!  Calling getElt on an unconverted array is not thread safe
3147        int error=0;
3148        if ((getDataPointRank()!=0))
3149        {
3150            throw DataException("Input to 1D interpolation must be scalar");
3151        }
3152        if (table.getRank()!=1)
3153        {
3154            throw DataException("Table for 1D interpolation must be 1D");
3155        }
3156        if (Astep<=0)
3157        {
3158            throw DataException("Astep must be positive");
3159        }
3160        if (!isExpanded())
3161        {
3162            expand();
3163        }
3164        Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);
3165        try
3166        {
3167            int numpts=getNumDataPoints();
3168            const DataVector& adat=getReady()->getVectorRO();
3169            DataVector& rdat=res.getReady()->getVectorRW();
3170            int twidth=table.getShape()[0]-1;
3171            bool haserror=false;
3172            int l=0;
3173            #pragma omp parallel for private(l) schedule(static)
3174            for (l=0;l<numpts; ++l)
3175            {
3176              #pragma omp flush(haserror)       // In case haserror was in register
3177              if (!haserror)        
3178              {
3179               int lerror=0;
3180               try
3181               {
3182                 do         // so we can use break
3183                 {
3184                double a=adat[l];
3185                int x=static_cast<int>(((a-Amin)/Astep));
3186                if (check_boundaries) {
3187                    if ((a<Amin) || (x<0))
3188                    {
3189                        lerror=1;
3190                        break;  
3191                    }
3192                    if (a>Amin+Astep*twidth)
3193                    {
3194                            lerror=4;
3195                        break;
3196                    }
3197                }
3198                if (x<0) x=0;
3199                if (x>twidth) x=twidth;
3200    
3201                if (x==twidth)          // value is on the far end of the table
3202                {
3203                    double e=table.getElt(x);
3204                    if (e>undef)
3205                    {
3206                        lerror=2;
3207                        break;
3208                    }
3209                    rdat[l]=e;
3210                }
3211                else        // x and y are in bounds
3212                {
3213                    double e=table.getElt(x);
3214                    double w=table.getElt(x+1);
3215                    if ((e>undef) || (w>undef))
3216                    {
3217                        lerror=2;
3218                        break;
3219                    }
3220            // map x*Astep <= a << (x+1)*Astep to [-1,1]
3221                    double la = 2.0*(a-Amin-(x*Astep))/Astep-1;
3222                    rdat[l]=((1-la)*e + (1+la)*w)/2;
3223                }  
3224                  } while (false);
3225                } catch (DataException d)
3226                {
3227                    lerror=3;
3228                }
3229                if (lerror!=0)
3230                {
3231                #pragma omp critical    // Doco says there is a flush associated with critical
3232                {
3233                    haserror=true;  // We only care that one error is recorded. We don't care which
3234                    error=lerror;   // one
3235                }
3236                }
3237              } // if (!error)
3238            }   // parallelised for
3239        } catch (DataException d)
3240        {
3241            error=3;        // this is outside the parallel region so assign directly
3242        }
3243    #ifdef PASO_MPI
3244        int rerror=0;
3245        MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );
3246        error=rerror;
3247    #endif
3248        if (error)
3249        {
3250            switch (error)
3251            {
3252                case 1: throw DataException("Value below lower table range.");
3253                case 2: throw DataException("Interpolated value too large");
3254                case 4: throw DataException("Value greater than upper table range.");
3255                default:
3256                    throw DataException("Unknown error in interpolation");      
3257            }
3258        }
3259        return res;
3260    }
3261    
3262            
3263    Data
3264    Data::interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
3265                           double undef, Data& B, double Bmin, double Bstep, bool check_boundaries)
3266    {
3267        table.convertArray();       // critical!   Calling getElt on an unconverted array is not thread safe
3268        int error=0;
3269        if ((getDataPointRank()!=0) || (B.getDataPointRank()!=0))
3270        {
3271            throw DataException("Inputs to 2D interpolation must be scalar");
3272        }
3273        if (table.getRank()!=2)
3274        {
3275        throw DataException("Table for 2D interpolation must be 2D");
3276        }
3277        if ((Astep<=0) || (Bstep<=0))
3278        {
3279        throw DataException("Astep and Bstep must be postive");
3280        }
3281        if (getFunctionSpace()!=B.getFunctionSpace())
3282        {
3283        Data n=B.interpolate(getFunctionSpace());
3284        return interpolateFromTable2D(table, Amin, Astep, undef,
3285            n , Bmin, Bstep, check_boundaries);
3286        }
3287        if (!isExpanded())
3288        {
3289        expand();
3290        }
3291        if (!B.isExpanded())
3292        {
3293        B.expand();
3294        }
3295        Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);
3296        try
3297        {
3298        int numpts=getNumDataPoints();
3299        const DataVector& adat=getReady()->getVectorRO();
3300        const DataVector& bdat=B.getReady()->getVectorRO();
3301        DataVector& rdat=res.getReady()->getVectorRW();
3302        const DataTypes::ShapeType& ts=table.getShape();
3303        int twx=ts[0]-1;    // table width x
3304        int twy=ts[1]-1;    // table width y
3305        bool haserror=false;
3306        int l=0;
3307        #pragma omp parallel for private(l) schedule(static)
3308        for (l=0; l<numpts; ++l)
3309        {
3310           #pragma omp flush(haserror)      // In case haserror was in register
3311           if (!haserror)      
3312           {
3313             int lerror=0;
3314             try
3315             {
3316               do
3317               {
3318            double a=adat[l];
3319            double b=bdat[l];
3320            int x=static_cast<int>(((a-Amin)/Astep));
3321            int y=static_cast<int>(((b-Bmin)/Bstep));
3322                    if (check_boundaries) {
3323                    if ( (a<Amin) || (b<Bmin) || (x<0) || (y<0) )
3324                    {
3325                    #pragma omp critical
3326                    {
3327                        lerror=1;
3328                    }
3329                    break;  
3330                    }
3331                    if ( (a>Amin+Astep*twx) || (b>Bmin+Bstep*twy) )
3332                    {
3333                    #pragma omp critical
3334                    {
3335                        lerror=4;
3336                    }
3337                    break;
3338                    }
3339                    }
3340                    if (x<0) x=0;
3341                    if (y<0) y=0;
3342                    if (x>twx) x=twx;
3343                    if (y>twx) y=twy;
3344    
3345                    if (x == twx ) {
3346                         if (y == twy ) {
3347                     double sw=table.getElt(x,y);
3348                     if ((sw>undef))
3349                     {
3350                     #pragma omp critical
3351                     {
3352                         lerror=2;
3353                     }
3354                     break;
3355                     }
3356                     rdat[l]=sw;
3357    
3358                         } else {
3359                     double sw=table.getElt(x,y);
3360                     double nw=table.getElt(x,y+1);
3361                     if ((sw>undef) || (nw>undef))
3362                     {
3363                     #pragma omp critical
3364                     {
3365                        lerror=2;
3366                     }
3367                     break;
3368                     }
3369                     double lb = 2.0*(b-Bmin-(y*Bstep))/Bstep-1;
3370                     rdat[l]=((1-lb)*sw + (1+lb)*nw )/2.;
3371    
3372                         }
3373                    } else {
3374                         if (y == twy ) {
3375                     double sw=table.getElt(x,y);
3376                     double se=table.getElt(x+1,y);
3377                     if ((sw>undef) || (se>undef) )
3378                     {
3379                     #pragma omp critical
3380                     {
3381                        lerror=2;
3382                     }
3383                     break;
3384                     }
3385                     double la = 2.0*(a-Amin-(x*Astep))/Astep-1;
3386                     rdat[l]=((1-la)*sw + (1+la)*se )/2;
3387    
3388                         } else {
3389                     double sw=table.getElt(x,y);
3390                     double nw=table.getElt(x,y+1);
3391                     double se=table.getElt(x+1,y);
3392                     double ne=table.getElt(x+1,y+1);
3393                     if ((sw>undef) || (nw>undef) || (se>undef) || (ne>undef))
3394                     {
3395                    #pragma omp critical
3396                    {
3397                         lerror=2;
3398                    }
3399                    break;
3400                     }
3401                     // map x*Astep <= a << (x+1)*Astep to [-1,1]
3402                     // same with b
3403                     double la = 2.0*(a-Amin-(x*Astep))/Astep-1;
3404                     double lb = 2.0*(b-Bmin-(y*Bstep))/Bstep-1;
3405                     rdat[l]=((1-la)*(1-lb)*sw + (1-la)*(1+lb)*nw +
3406                          (1+la)*(1-lb)*se + (1+la)*(1+lb)*ne)/4;
3407                         }
3408                    }
3409               } while (false);
3410              } catch (DataException d)
3411              {
3412            lerror=3;
3413              }
3414              if (lerror!=0)
3415              {
3416                #pragma omp critical    // Doco says there is a flush associated with critical
3417                {
3418                    haserror=true;  // We only care that one error is recorded. We don't care which
3419                    error=lerror;   // one
3420                }      
3421              }
3422            }  // if (!error)
3423        }   // parallel for
3424        } catch (DataException d)
3425        {
3426            error=3;
3427        }
3428    #ifdef PASO_MPI
3429        int rerror=0;
3430        MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );
3431        error=rerror;
3432    #endif
3433        if (error)
3434        {
3435        switch (error)
3436        {
3437                case 1: throw DataException("Value below lower table range.");
3438                case 2: throw DataException("Interpolated value too large");
3439                case 4: throw DataException("Value greater than upper table range.");
3440                        default:
3441                          throw DataException("Unknown error in interpolation");        
3442        }
3443        }
3444        return res;
3445  }  }
3446    
3447    
# Line 3000  Data::print() Line 3457  Data::print()
3457    {    {
3458      printf( "[%6d]", i );      printf( "[%6d]", i );
3459      for( j=0; j<getNumDataPointsPerSample(); j++ )      for( j=0; j<getNumDataPointsPerSample(); j++ )
3460        printf( "\t%10.7g", (getSampleData(i))[j] );        printf( "\t%10.7g", (getSampleDataRW(i))[j] );    // doesn't really need RW access
3461      printf( "\n" );      printf( "\n" );
3462    }    }
3463  }  }
# Line 3020  Data::dump(const std::string fileName) c Line 3477  Data::dump(const std::string fileName) c
3477            return m_data->dump(fileName);            return m_data->dump(fileName);
3478      }      }
3479    }    }
3480    catch (exception& e)    catch (std::exception& e)
3481    {    {
3482          cout << e.what() << endl;          cout << e.what() << endl;
3483    }    }
# Line 3062  Data::get_MPIComm() const Line 3519  Data::get_MPIComm() const
3519  #endif  #endif
3520  }  }
3521    
   

Legend:
Removed from v.2049  
changed lines
  Added in v.2735

  ViewVC Help
Powered by ViewVC 1.1.26