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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2635 - (hide annotations)
Thu Aug 27 04:54:41 2009 UTC (10 years, 4 months ago) by jfenwick
File size: 89284 byte(s)
A bunch of changes related to saveDataCSV.
[Not completed or unit tested yet]

Added saveDataCSV to util.py
AbstractDomain (and MeshAdapter) have a commonFunctionSpace method to 
take a group of FunctionSpaces and return something they can all be interpolated to.

Added pointToStream() in DataTypes to help print points.

added actsConstant() to data - required because DataConstant doesn't store samples the same way other Data do.
1 jgs 480
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1811 * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 ksteube 1312
14 ksteube 1811
15 jgs 474 #include "Data.h"
16 jgs 94
17 jgs 480 #include "DataExpanded.h"
18     #include "DataConstant.h"
19     #include "DataTagged.h"
20     #include "DataEmpty.h"
21 jfenwick 2005 #include "DataLazy.h"
22 jgs 480 #include "FunctionSpaceFactory.h"
23     #include "AbstractContinuousDomain.h"
24     #include "UnaryFuncs.h"
25 jfenwick 1796 #include "FunctionSpaceException.h"
26 jfenwick 1897 #include "EscriptParams.h"
27 jfenwick 1796
28 ksteube 1312 extern "C" {
29 phornby 2078 #include "esysUtils/blocktimer.h"
30 ksteube 1312 }
31 jgs 480
32 jgs 119 #include <fstream>
33 jgs 94 #include <algorithm>
34     #include <vector>
35     #include <functional>
36 jfenwick 2086 #include <sstream> // so we can throw messages about ranks
37 jgs 94
38 jgs 153 #include <boost/python/dict.hpp>
39 jgs 94 #include <boost/python/extract.hpp>
40     #include <boost/python/long.hpp>
41 jfenwick 2271 #include "WrappedArray.h"
42 jgs 94
43     using namespace std;
44     using namespace boost::python;
45     using namespace boost;
46     using namespace escript;
47    
48 jfenwick 2005 // ensure the current object is not a DataLazy
49     // The idea was that we could add an optional warning whenever a resolve is forced
50 jfenwick 2271 // #define forceResolve() if (isLazy()) {#resolve();}
51    
52 jfenwick 2146 #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 jfenwick 2005
64 jfenwick 2496 #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 jfenwick 2146 #define MAKELAZYBINSELF(R,X) if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
71     {\
72     DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
73 jfenwick 2271 /* m_data=c->getPtr();*/ set_m_data(c->getPtr());\
74 jfenwick 2146 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 jfenwick 2271 // Do not use the following unless you want to make copies on assignment rather than
91     // share data. CopyOnWrite should make this unnescessary.
92     // #define ASSIGNMENT_MEANS_DEEPCOPY
93    
94     namespace
95     {
96    
97 jfenwick 2459 template <class ARR>
98     inline
99     boost::python::tuple
100     pointToTuple1(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
101     {
102     using namespace boost::python;
103     using boost::python::tuple;
104     using boost::python::list;
105    
106     list l;
107     unsigned int dim=shape[0];
108     for (size_t i=0;i<dim;++i)
109     {
110     l.append(v[i+offset]);
111     }
112     return tuple(l);
113     }
114    
115     template <class ARR>
116     inline
117     boost::python::tuple
118     pointToTuple2(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
119     {
120     using namespace boost::python;
121     using boost::python::tuple;
122     using boost::python::list;
123    
124     unsigned int shape0=shape[0];
125     unsigned int shape1=shape[1];
126     list lj;
127     for (size_t j=0;j<shape0;++j)
128     {
129     list li;
130     for (size_t i=0;i<shape1;++i)
131     {
132     li.append(v[offset+DataTypes::getRelIndex(shape,j,i)]);
133     }
134     lj.append(tuple(li));
135     }
136     return tuple(lj);
137     }
138    
139     template <class ARR>
140     inline
141     boost::python::tuple
142     pointToTuple3(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
143     {
144     using namespace boost::python;
145     using boost::python::tuple;
146     using boost::python::list;
147    
148     unsigned int shape0=shape[0];
149     unsigned int shape1=shape[1];
150     unsigned int shape2=shape[2];
151    
152     list lk;
153     for (size_t k=0;k<shape0;++k)
154     {
155     list lj;
156     for (size_t j=0;j<shape1;++j)
157     {
158     list li;
159     for (size_t i=0;i<shape2;++i)
160     {
161 jfenwick 2482 li.append(v[offset+DataTypes::getRelIndex(shape,k,j,i)]);
162 jfenwick 2459 }
163     lj.append(tuple(li));
164     }
165     lk.append(tuple(lj));
166     }
167     return tuple(lk);
168     }
169    
170     template <class ARR>
171     inline
172     boost::python::tuple
173     pointToTuple4(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
174     {
175     using namespace boost::python;
176     using boost::python::tuple;
177     using boost::python::list;
178    
179     unsigned int shape0=shape[0];
180     unsigned int shape1=shape[1];
181     unsigned int shape2=shape[2];
182     unsigned int shape3=shape[3];
183    
184     list ll;
185     for (size_t l=0;l<shape0;++l)
186     {
187     list lk;
188     for (size_t k=0;k<shape1;++k)
189     {
190     list lj;
191     for (size_t j=0;j<shape2;++j)
192     {
193     list li;
194     for (size_t i=0;i<shape3;++i)
195     {
196 jfenwick 2482 li.append(v[offset+DataTypes::getRelIndex(shape,l,k,j,i)]);
197 jfenwick 2459 }
198     lj.append(tuple(li));
199     }
200     lk.append(tuple(lj));
201     }
202     ll.append(tuple(lk));
203     }
204     return tuple(ll);
205     }
206    
207    
208 jfenwick 2271 // This should be safer once the DataC RO changes have been brought in
209     template <class ARR>
210     boost::python::tuple
211     pointToTuple( const DataTypes::ShapeType& shape,ARR v)
212     {
213     using namespace boost::python;
214     using boost::python::list;
215     int rank=shape.size();
216     if (rank==0)
217     {
218     return make_tuple(v[0]);
219     }
220     else if (rank==1)
221     {
222 jfenwick 2459 return pointToTuple1(shape,v,0);
223 jfenwick 2271 }
224     else if (rank==2)
225     {
226 jfenwick 2459 return pointToTuple2(shape,v,0);
227 jfenwick 2271 }
228     else if (rank==3)
229     {
230 jfenwick 2459 return pointToTuple3(shape,v,0);
231 jfenwick 2271 }
232     else if (rank==4)
233     {
234 jfenwick 2459 return pointToTuple4(shape,v,0);
235 jfenwick 2271 }
236     else
237     throw DataException("Unknown rank in pointToTuple.");
238     }
239    
240     } // anonymous namespace
241    
242 jgs 94 Data::Data()
243 jfenwick 2271 : m_shared(false), m_lazy(false)
244 jgs 94 {
245     //
246     // Default data is type DataEmpty
247     DataAbstract* temp=new DataEmpty();
248 jfenwick 2271 // m_data=temp->getPtr();
249     set_m_data(temp->getPtr());
250 gross 783 m_protected=false;
251 jgs 94 }
252    
253     Data::Data(double value,
254     const tuple& shape,
255     const FunctionSpace& what,
256     bool expanded)
257 jfenwick 2271 : m_shared(false), m_lazy(false)
258 jgs 94 {
259 jfenwick 1796 DataTypes::ShapeType dataPointShape;
260 jgs 94 for (int i = 0; i < shape.attr("__len__")(); ++i) {
261     dataPointShape.push_back(extract<const int>(shape[i]));
262     }
263 matt 1319
264 jfenwick 1796 int len = DataTypes::noValues(dataPointShape);
265 matt 1319 DataVector temp_data(len,value,len);
266 jfenwick 1796 initialise(temp_data, dataPointShape, what, expanded);
267 gross 783 m_protected=false;
268 jgs 94 }
269    
270     Data::Data(double value,
271 jfenwick 1796 const DataTypes::ShapeType& dataPointShape,
272 jgs 94 const FunctionSpace& what,
273     bool expanded)
274 jfenwick 2271 : m_shared(false), m_lazy(false)
275 jgs 94 {
276 jfenwick 1796 int len = DataTypes::noValues(dataPointShape);
277 matt 1319 DataVector temp_data(len,value,len);
278 jfenwick 1796 initialise(temp_data, dataPointShape, what, expanded);
279 gross 783 m_protected=false;
280 jgs 94 }
281    
282 jgs 102 Data::Data(const Data& inData)
283 jfenwick 2271 : m_shared(false), m_lazy(false)
284 jgs 94 {
285 jfenwick 2271 set_m_data(inData.m_data);
286 gross 783 m_protected=inData.isProtected();
287 jgs 94 }
288    
289 jfenwick 1796
290 jgs 94 Data::Data(const Data& inData,
291 jfenwick 1796 const DataTypes::RegionType& region)
292 jfenwick 2271 : m_shared(false), m_lazy(false)
293 jgs 94 {
294 jfenwick 2005 DataAbstract_ptr dat=inData.m_data;
295     if (inData.isLazy())
296     {
297     dat=inData.m_data->resolve();
298     }
299     else
300     {
301     dat=inData.m_data;
302     }
303 jgs 94 //
304 jgs 102 // Create Data which is a slice of another Data
305 jfenwick 2005 DataAbstract* tmp = dat->getSlice(region);
306 jfenwick 2271 set_m_data(DataAbstract_ptr(tmp));
307 gross 783 m_protected=false;
308 jfenwick 2271
309 jgs 94 }
310    
311     Data::Data(const Data& inData,
312     const FunctionSpace& functionspace)
313 jfenwick 2271 : m_shared(false), m_lazy(false)
314 jgs 94 {
315 jfenwick 1803 if (inData.isEmpty())
316     {
317     throw DataException("Error - will not interpolate for instances of DataEmpty.");
318     }
319 jgs 94 if (inData.getFunctionSpace()==functionspace) {
320 jfenwick 2271 set_m_data(inData.m_data);
321 jfenwick 2005 }
322     else
323     {
324    
325     if (inData.isConstant()) { // for a constant function, we just need to use the new function space
326     if (!inData.probeInterpolation(functionspace))
327     { // Even though this is constant, we still need to check whether interpolation is allowed
328 jfenwick 2087 throw FunctionSpaceException("Cannot interpolate across to the domain of the specified FunctionSpace. (DataConstant)");
329 jfenwick 2005 }
330     // if the data is not lazy, this will just be a cast to DataReady
331     DataReady_ptr dr=inData.m_data->resolve();
332 jfenwick 2271 DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVectorRO());
333     // m_data=DataAbstract_ptr(dc);
334     set_m_data(DataAbstract_ptr(dc));
335 jgs 94 } else {
336 jfenwick 2005 Data tmp(0,inData.getDataPointShape(),functionspace,true);
337     // Note: Must use a reference or pointer to a derived object
338     // in order to get polymorphic behaviour. Shouldn't really
339     // be able to create an instance of AbstractDomain but that was done
340     // as a boost:python work around which may no longer be required.
341     /*const AbstractDomain& inDataDomain=inData.getDomain();*/
342     const_Domain_ptr inDataDomain=inData.getDomain();
343     if (inDataDomain==functionspace.getDomain()) {
344     inDataDomain->interpolateOnDomain(tmp,inData);
345     } else {
346     inDataDomain->interpolateACross(tmp,inData);
347     }
348 jfenwick 2271 // m_data=tmp.m_data;
349     set_m_data(tmp.m_data);
350 jgs 94 }
351     }
352 gross 783 m_protected=false;
353 jgs 94 }
354    
355 jfenwick 1796 Data::Data(DataAbstract* underlyingdata)
356 jfenwick 2271 : m_shared(false), m_lazy(false)
357 jgs 94 {
358 jfenwick 2271 set_m_data(underlyingdata->getPtr());
359 jfenwick 1796 m_protected=false;
360 jgs 94 }
361    
362 jfenwick 2005 Data::Data(DataAbstract_ptr underlyingdata)
363 jfenwick 2271 : m_shared(false), m_lazy(false)
364 jfenwick 2005 {
365 jfenwick 2271 set_m_data(underlyingdata);
366 jfenwick 2005 m_protected=false;
367     }
368    
369 jfenwick 1796 Data::Data(const DataTypes::ValueType& value,
370     const DataTypes::ShapeType& shape,
371     const FunctionSpace& what,
372     bool expanded)
373 jfenwick 2271 : m_shared(false), m_lazy(false)
374 jfenwick 1796 {
375     initialise(value,shape,what,expanded);
376     m_protected=false;
377 jgs 94 }
378    
379 jfenwick 1796
380 jgs 94 Data::Data(const object& value,
381     const FunctionSpace& what,
382     bool expanded)
383 jfenwick 2271 : m_shared(false), m_lazy(false)
384 jgs 94 {
385 jfenwick 2271 WrappedArray w(value);
386     initialise(w,what,expanded);
387 gross 783 m_protected=false;
388 jgs 94 }
389    
390 matt 1319
391 jgs 94 Data::Data(const object& value,
392     const Data& other)
393 jfenwick 2271 : m_shared(false), m_lazy(false)
394 jgs 94 {
395 jfenwick 2271 WrappedArray w(value);
396 matt 1319
397 jfenwick 2458 // extract the shape of the array
398 jfenwick 2271 const DataTypes::ShapeType& tempShape=w.getShape();
399     if (w.getRank()==0) {
400 matt 1319
401    
402 jfenwick 1796 // get the space for the data vector
403     int len1 = DataTypes::noValues(tempShape);
404     DataVector temp_data(len1, 0.0, len1);
405 jfenwick 2271 temp_data.copyFromArray(w,1);
406 jfenwick 1796
407     int len = DataTypes::noValues(other.getDataPointShape());
408    
409 jfenwick 2271 DataVector temp2_data(len, temp_data[0], len);
410 jfenwick 1796 DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
411 jfenwick 2271 // m_data=DataAbstract_ptr(t);
412     set_m_data(DataAbstract_ptr(t));
413 jfenwick 1796
414 jgs 94 } else {
415     //
416     // Create a DataConstant with the same sample shape as other
417 jfenwick 2271 DataConstant* t=new DataConstant(w,other.getFunctionSpace());
418     // m_data=DataAbstract_ptr(t);
419     set_m_data(DataAbstract_ptr(t));
420 jgs 94 }
421 gross 783 m_protected=false;
422 jgs 94 }
423    
424 jgs 151 Data::~Data()
425     {
426 jfenwick 2271 set_m_data(DataAbstract_ptr());
427 jgs 151 }
428    
429 jfenwick 1796
430 jfenwick 2271 // only call in thread safe contexts.
431     // This method should be atomic
432     void Data::set_m_data(DataAbstract_ptr p)
433     {
434     if (m_data.get()!=0) // release old ownership
435     {
436     m_data->removeOwner(this);
437     }
438     if (p.get()!=0)
439     {
440     m_data=p;
441     m_data->addOwner(this);
442     m_shared=m_data->isShared();
443     m_lazy=m_data->isLazy();
444     }
445     }
446 jfenwick 1796
447 jfenwick 2271 void Data::initialise(const WrappedArray& value,
448 jfenwick 1796 const FunctionSpace& what,
449     bool expanded)
450     {
451     //
452     // Construct a Data object of the appropriate type.
453     // Construct the object first as there seems to be a bug which causes
454     // undefined behaviour if an exception is thrown during construction
455     // within the shared_ptr constructor.
456     if (expanded) {
457     DataAbstract* temp=new DataExpanded(value, what);
458 jfenwick 2271 // m_data=temp->getPtr();
459     set_m_data(temp->getPtr());
460 jfenwick 1796 } else {
461     DataAbstract* temp=new DataConstant(value, what);
462 jfenwick 2271 // m_data=temp->getPtr();
463     set_m_data(temp->getPtr());
464 jfenwick 1796 }
465     }
466    
467    
468     void
469     Data::initialise(const DataTypes::ValueType& value,
470     const DataTypes::ShapeType& shape,
471     const FunctionSpace& what,
472     bool expanded)
473     {
474     //
475     // Construct a Data object of the appropriate type.
476     // Construct the object first as there seems to be a bug which causes
477     // undefined behaviour if an exception is thrown during construction
478     // within the shared_ptr constructor.
479     if (expanded) {
480     DataAbstract* temp=new DataExpanded(what, shape, value);
481 jfenwick 2271 // m_data=temp->getPtr();
482     set_m_data(temp->getPtr());
483 jfenwick 1796 } else {
484     DataAbstract* temp=new DataConstant(what, shape, value);
485 jfenwick 2271 // m_data=temp->getPtr();
486     set_m_data(temp->getPtr());
487 jfenwick 1796 }
488     }
489    
490    
491 jgs 94 escriptDataC
492     Data::getDataC()
493     {
494     escriptDataC temp;
495     temp.m_dataPtr=(void*)this;
496     return temp;
497     }
498    
499     escriptDataC
500     Data::getDataC() const
501     {
502     escriptDataC temp;
503     temp.m_dataPtr=(void*)this;
504     return temp;
505     }
506    
507 jfenwick 2271 size_t
508     Data::getSampleBufferSize() const
509     {
510     return m_data->getSampleBufferSize();
511     }
512    
513    
514 jgs 121 const boost::python::tuple
515 jgs 94 Data::getShapeTuple() const
516     {
517 jfenwick 1796 const DataTypes::ShapeType& shape=getDataPointShape();
518 jgs 94 switch(getDataPointRank()) {
519     case 0:
520     return make_tuple();
521     case 1:
522     return make_tuple(long_(shape[0]));
523     case 2:
524     return make_tuple(long_(shape[0]),long_(shape[1]));
525     case 3:
526     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]));
527     case 4:
528     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]),long_(shape[3]));
529     default:
530     throw DataException("Error - illegal Data rank.");
531     }
532     }
533 jfenwick 1799
534    
535     // The different name is needed because boost has trouble with overloaded functions.
536     // It can't work out what type the function is based soley on its name.
537     // There are ways to fix this involving creating function pointer variables for each form
538     // but there doesn't seem to be a need given that the methods have the same name from the python point of view
539 jfenwick 2105 Data
540 jfenwick 1799 Data::copySelf()
541     {
542     DataAbstract* temp=m_data->deepCopy();
543 jfenwick 2105 return Data(temp);
544 jfenwick 1799 }
545    
546 jgs 94 void
547     Data::copy(const Data& other)
548     {
549 jfenwick 1799 DataAbstract* temp=other.m_data->deepCopy();
550 jfenwick 1872 DataAbstract_ptr p=temp->getPtr();
551 jfenwick 2271 // m_data=p;
552     set_m_data(p);
553 jgs 94 }
554    
555 gross 1118
556 jfenwick 2005 Data
557     Data::delay()
558     {
559     DataLazy* dl=new DataLazy(m_data);
560     return Data(dl);
561     }
562    
563 jgs 94 void
564 jfenwick 2005 Data::delaySelf()
565     {
566     if (!isLazy())
567     {
568 jfenwick 2271 // m_data=(new DataLazy(m_data))->getPtr();
569     set_m_data((new DataLazy(m_data))->getPtr());
570 jfenwick 2005 }
571     }
572    
573 jfenwick 2271
574     // For lazy data, it would seem that DataTagged will need to be treated differently since even after setting all tags
575     // to zero, all the tags from all the DataTags would be in the result.
576     // However since they all have the same value (0) whether they are there or not should not matter.
577     // So I have decided that for all types this method will create a constant 0.
578     // It can be promoted up as required.
579     // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
580     // but we can deal with that if it arises.
581     //
582 jfenwick 2005 void
583 gross 1118 Data::setToZero()
584     {
585 jfenwick 1803 if (isEmpty())
586 gross 1118 {
587 jfenwick 1803 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
588     }
589 jfenwick 2271 if (isLazy())
590     {
591     DataTypes::ValueType v(getNoValues(),0);
592     DataConstant* dc=new DataConstant(getFunctionSpace(),getDataPointShape(),v);
593     DataLazy* dl=new DataLazy(dc->getPtr());
594     set_m_data(dl->getPtr());
595     }
596     else
597     {
598     exclusiveWrite();
599     m_data->setToZero();
600     }
601 gross 1118 }
602    
603 jfenwick 2271
604 gross 1118 void
605 jgs 94 Data::copyWithMask(const Data& other,
606     const Data& mask)
607     {
608 jfenwick 1855 // 1. Interpolate if required so all Datas use the same FS as this
609     // 2. Tag or Expand so that all Data's are the same type
610     // 3. Iterate over the data vectors copying values where mask is >0
611 jfenwick 1803 if (other.isEmpty() || mask.isEmpty())
612     {
613     throw DataException("Error - copyWithMask not permitted using instances of DataEmpty.");
614     }
615 jfenwick 1855 Data other2(other);
616     Data mask2(mask);
617 jfenwick 2005 other2.resolve();
618     mask2.resolve();
619     this->resolve();
620 jfenwick 1855 FunctionSpace myFS=getFunctionSpace();
621     FunctionSpace oFS=other2.getFunctionSpace();
622     FunctionSpace mFS=mask2.getFunctionSpace();
623     if (oFS!=myFS)
624     {
625     if (other2.probeInterpolation(myFS))
626     {
627     other2=other2.interpolate(myFS);
628     }
629     else
630     {
631     throw DataException("Error - copyWithMask: other FunctionSpace is not compatible with this one.");
632     }
633     }
634     if (mFS!=myFS)
635     {
636     if (mask2.probeInterpolation(myFS))
637     {
638     mask2=mask2.interpolate(myFS);
639     }
640     else
641     {
642     throw DataException("Error - copyWithMask: mask FunctionSpace is not compatible with this one.");
643     }
644     }
645     // Ensure that all args have the same type
646     if (this->isExpanded() || mask2.isExpanded() || other2.isExpanded())
647     {
648     this->expand();
649     other2.expand();
650     mask2.expand();
651     }
652     else if (this->isTagged() || mask2.isTagged() || other2.isTagged())
653     {
654     this->tag();
655     other2.tag();
656     mask2.tag();
657     }
658     else if (this->isConstant() && mask2.isConstant() && other2.isConstant())
659     {
660     }
661     else
662     {
663     throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
664     }
665 jfenwick 2246 unsigned int selfrank=getDataPointRank();
666     unsigned int otherrank=other2.getDataPointRank();
667     unsigned int maskrank=mask2.getDataPointRank();
668     if ((selfrank==0) && (otherrank>0 || maskrank>0))
669     {
670     // to get here we must be copying from a large object into a scalar
671     // I am not allowing this.
672     // If you are calling copyWithMask then you are considering keeping some existing values
673     // and so I'm going to assume that you don't want your data objects getting a new shape.
674     throw DataException("Attempt to copyWithMask into a scalar from an object or mask with rank>0.");
675     }
676 jfenwick 2271 exclusiveWrite();
677 jfenwick 1855 // Now we iterate over the elements
678 jfenwick 2271 DataVector& self=getReady()->getVectorRW();;
679     const DataVector& ovec=other2.getReadyPtr()->getVectorRO();
680     const DataVector& mvec=mask2.getReadyPtr()->getVectorRO();
681 jfenwick 2246
682     if ((selfrank>0) && (otherrank==0) &&(maskrank==0))
683 jfenwick 1855 {
684 jfenwick 2246 // Not allowing this combination.
685     // it is not clear what the rank of the target should be.
686     // Should it be filled with the scalar (rank stays the same);
687     // or should the target object be reshaped to be a scalar as well.
688     throw DataException("Attempt to copyWithMask from scalar mask and data into non-scalar target.");
689 jfenwick 1855 }
690 jfenwick 2246 if ((selfrank>0) && (otherrank>0) &&(maskrank==0))
691     {
692     if (mvec[0]>0) // copy whole object if scalar is >0
693     {
694     copy(other);
695     }
696     return;
697     }
698     if (isTagged()) // so all objects involved will also be tagged
699     {
700     // note the !
701     if (!((getDataPointShape()==mask2.getDataPointShape()) &&
702     ((other2.getDataPointShape()==mask2.getDataPointShape()) || (otherrank==0))))
703     {
704     throw DataException("copyWithMask, shape mismatch.");
705     }
706    
707     // We need to consider the possibility that tags are missing or in the wrong order
708     // My guiding assumption here is: All tagged Datas are assumed to have the default value for
709     // all tags which are not explicitly defined
710    
711     const DataTagged* mptr=dynamic_cast<const DataTagged*>(mask2.m_data.get());
712     const DataTagged* optr=dynamic_cast<const DataTagged*>(other2.m_data.get());
713     DataTagged* tptr=dynamic_cast<DataTagged*>(m_data.get());
714    
715     // first, add any tags missing from other or mask
716     const DataTagged::DataMapType& olookup=optr->getTagLookup();
717     const DataTagged::DataMapType& mlookup=mptr->getTagLookup();
718 jfenwick 2260 const DataTagged::DataMapType& tlookup=tptr->getTagLookup();
719 jfenwick 2246 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
720     for (i=olookup.begin();i!=olookup.end();i++)
721     {
722     tptr->addTag(i->first);
723     }
724     for (i=mlookup.begin();i!=mlookup.end();i++) {
725     tptr->addTag(i->first);
726     }
727     // now we know that *this has all the required tags but they aren't guaranteed to be in
728     // the same order
729    
730     // There are two possibilities: 1. all objects have the same rank. 2. other is a scalar
731     if ((selfrank==otherrank) && (otherrank==maskrank))
732     {
733 jfenwick 2260 for (i=tlookup.begin();i!=tlookup.end();i++)
734 jfenwick 2246 {
735     // get the target offset
736     DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
737     DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
738     DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
739 jfenwick 2260 for (int j=0;j<getDataPointSize();++j)
740 jfenwick 2246 {
741 jfenwick 2260 if (mvec[j+moff]>0)
742 jfenwick 2246 {
743 jfenwick 2260 self[j+toff]=ovec[j+ooff];
744 jfenwick 2246 }
745     }
746     }
747 jfenwick 2260 // now for the default value
748     for (int j=0;j<getDataPointSize();++j)
749     {
750     if (mvec[j+mptr->getDefaultOffset()]>0)
751     {
752     self[j+tptr->getDefaultOffset()]=ovec[j+optr->getDefaultOffset()];
753     }
754     }
755 jfenwick 2246 }
756     else // other is a scalar
757     {
758 jfenwick 2260 for (i=tlookup.begin();i!=tlookup.end();i++)
759 jfenwick 2246 {
760     // get the target offset
761     DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
762     DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
763     DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
764 jfenwick 2260 for (int j=0;j<getDataPointSize();++j)
765 jfenwick 2246 {
766 jfenwick 2260 if (mvec[j+moff]>0)
767 jfenwick 2246 {
768 jfenwick 2260 self[j+toff]=ovec[ooff];
769 jfenwick 2246 }
770     }
771     }
772 jfenwick 2260 // now for the default value
773     for (int j=0;j<getDataPointSize();++j)
774     {
775     if (mvec[j+mptr->getDefaultOffset()]>0)
776     {
777     self[j+tptr->getDefaultOffset()]=ovec[0];
778     }
779     }
780 jfenwick 2246 }
781    
782     return; // ugly
783     }
784     // mixed scalar and non-scalar operation
785     if ((selfrank>0) && (otherrank==0) && (mask2.getDataPointShape()==getDataPointShape()))
786     {
787     size_t num_points=self.size();
788     // OPENMP 3.0 allows unsigned loop vars.
789     #if defined(_OPENMP) && (_OPENMP < 200805)
790     long i;
791     #else
792     size_t i;
793 jfenwick 2260 #endif
794     size_t psize=getDataPointSize();
795 jfenwick 2246 #pragma omp parallel for private(i) schedule(static)
796     for (i=0;i<num_points;++i)
797     {
798     if (mvec[i]>0)
799     {
800 jfenwick 2260 self[i]=ovec[i/psize]; // since this is expanded there is one scalar
801     } // dest point
802 jfenwick 2246 }
803     return; // ugly!
804     }
805     // tagged data is already taken care of so we only need to worry about shapes
806     // special cases with scalars are already dealt with so all we need to worry about is shape
807     if ((getDataPointShape()!=other2.getDataPointShape()) || getDataPointShape()!=mask2.getDataPointShape())
808     {
809     ostringstream oss;
810     oss <<"Error - size mismatch in arguments to copyWithMask.";
811     oss << "\nself_shape=" << DataTypes::shapeToString(getDataPointShape());
812     oss << " other2_shape=" << DataTypes::shapeToString(other2.getDataPointShape());
813     oss << " mask2_shape=" << DataTypes::shapeToString(mask2.getDataPointShape());
814     throw DataException(oss.str());
815     }
816 jfenwick 1855 size_t num_points=self.size();
817 phornby 1921
818     // OPENMP 3.0 allows unsigned loop vars.
819     #if defined(_OPENMP) && (_OPENMP < 200805)
820     long i;
821     #else
822 phornby 1918 size_t i;
823 phornby 1921 #endif
824 jfenwick 1855 #pragma omp parallel for private(i) schedule(static)
825     for (i=0;i<num_points;++i)
826     {
827     if (mvec[i]>0)
828     {
829     self[i]=ovec[i];
830     }
831     }
832     }
833 jgs 94
834     bool
835     Data::isExpanded() const
836     {
837     DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
838     return (temp!=0);
839     }
840    
841     bool
842 jfenwick 2271 Data::actsExpanded() const
843     {
844     return m_data->actsExpanded();
845    
846     }
847    
848    
849     bool
850 jgs 94 Data::isTagged() const
851     {
852     DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
853     return (temp!=0);
854     }
855    
856     bool
857     Data::isEmpty() const
858     {
859     DataEmpty* temp=dynamic_cast<DataEmpty*>(m_data.get());
860     return (temp!=0);
861     }
862    
863     bool
864     Data::isConstant() const
865     {
866     DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
867     return (temp!=0);
868     }
869    
870 jfenwick 2005 bool
871 jfenwick 2635 Data::actsConstant() const
872     {
873     return m_data->actsConstant();
874     }
875    
876    
877     bool
878 jfenwick 2005 Data::isLazy() const
879     {
880 jfenwick 2271 return m_lazy; // not asking m_data because we need to be able to ask this while m_data is changing
881 jfenwick 2005 }
882    
883     // at the moment this is synonymous with !isLazy() but that could change
884     bool
885     Data::isReady() const
886     {
887     return (dynamic_cast<DataReady*>(m_data.get())!=0);
888     }
889    
890    
891 jgs 94 void
892 ksteube 1312 Data::setProtection()
893     {
894 gross 783 m_protected=true;
895     }
896    
897     bool
898 ksteube 1312 Data::isProtected() const
899     {
900 gross 783 return m_protected;
901     }
902    
903    
904    
905     void
906 jgs 94 Data::expand()
907     {
908     if (isConstant()) {
909     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
910     DataAbstract* temp=new DataExpanded(*tempDataConst);
911 jfenwick 2271 // m_data=temp->getPtr();
912     set_m_data(temp->getPtr());
913 jgs 94 } else if (isTagged()) {
914     DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
915     DataAbstract* temp=new DataExpanded(*tempDataTag);
916 jfenwick 2271 // m_data=temp->getPtr();
917     set_m_data(temp->getPtr());
918 jgs 94 } else if (isExpanded()) {
919     //
920     // do nothing
921     } else if (isEmpty()) {
922     throw DataException("Error - Expansion of DataEmpty not possible.");
923 jfenwick 2005 } else if (isLazy()) {
924     resolve();
925     expand(); // resolve might not give us expanded data
926 jgs 94 } else {
927     throw DataException("Error - Expansion not implemented for this Data type.");
928     }
929     }
930    
931     void
932     Data::tag()
933     {
934     if (isConstant()) {
935     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
936     DataAbstract* temp=new DataTagged(*tempDataConst);
937 jfenwick 2271 // m_data=temp->getPtr();
938     set_m_data(temp->getPtr());
939 jgs 94 } else if (isTagged()) {
940     // do nothing
941     } else if (isExpanded()) {
942     throw DataException("Error - Creating tag data from DataExpanded not possible.");
943     } else if (isEmpty()) {
944     throw DataException("Error - Creating tag data from DataEmpty not possible.");
945 jfenwick 2005 } else if (isLazy()) {
946     DataAbstract_ptr res=m_data->resolve();
947     if (m_data->isExpanded())
948     {
949     throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");
950     }
951 jfenwick 2271 // m_data=res;
952     set_m_data(res);
953 jfenwick 2005 tag();
954 jgs 94 } else {
955     throw DataException("Error - Tagging not implemented for this Data type.");
956     }
957     }
958    
959 jfenwick 2005 void
960     Data::resolve()
961     {
962     if (isLazy())
963     {
964 jfenwick 2271 // m_data=m_data->resolve();
965     set_m_data(m_data->resolve());
966 jfenwick 2005 }
967     }
968    
969 jfenwick 2271 void
970     Data::requireWrite()
971     {
972     resolve();
973     exclusiveWrite();
974     }
975 jfenwick 2005
976 gross 854 Data
977     Data::oneOver() const
978 jgs 102 {
979 jfenwick 2146 MAKELAZYOP(RECIP)
980 matt 1334 return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
981 jgs 102 }
982    
983 jgs 94 Data
984 gross 698 Data::wherePositive() const
985 jgs 94 {
986 jfenwick 2146 MAKELAZYOP(GZ)
987 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
988 jgs 94 }
989    
990     Data
991 gross 698 Data::whereNegative() const
992 jgs 102 {
993 jfenwick 2146 MAKELAZYOP(LZ)
994 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
995 jgs 102 }
996    
997     Data
998 gross 698 Data::whereNonNegative() const
999 jgs 94 {
1000 jfenwick 2146 MAKELAZYOP(GEZ)
1001 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
1002 jgs 94 }
1003    
1004     Data
1005 gross 698 Data::whereNonPositive() const
1006 jgs 94 {
1007 jfenwick 2146 MAKELAZYOP(LEZ)
1008 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
1009 jgs 94 }
1010    
1011     Data
1012 jgs 571 Data::whereZero(double tol) const
1013 jgs 94 {
1014 jfenwick 2147 // Data dataAbs=abs();
1015     // return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
1016     MAKELAZYOPOFF(EZ,tol)
1017     return C_TensorUnaryOperation(*this, bind2nd(AbsLTE(),tol));
1018    
1019 jgs 94 }
1020    
1021     Data
1022 jgs 571 Data::whereNonZero(double tol) const
1023 jgs 102 {
1024 jfenwick 2147 // Data dataAbs=abs();
1025     // return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
1026     MAKELAZYOPOFF(NEZ,tol)
1027     return C_TensorUnaryOperation(*this, bind2nd(AbsGT(),tol));
1028    
1029 jgs 102 }
1030    
1031     Data
1032 jgs 94 Data::interpolate(const FunctionSpace& functionspace) const
1033     {
1034     return Data(*this,functionspace);
1035     }
1036    
1037     bool
1038     Data::probeInterpolation(const FunctionSpace& functionspace) const
1039     {
1040 jfenwick 2005 return getFunctionSpace().probeInterpolation(functionspace);
1041 jgs 94 }
1042    
1043     Data
1044     Data::gradOn(const FunctionSpace& functionspace) const
1045     {
1046 jfenwick 1803 if (isEmpty())
1047     {
1048     throw DataException("Error - operation not permitted on instances of DataEmpty.");
1049     }
1050 matt 1353 double blocktimer_start = blocktimer_time();
1051 jgs 94 if (functionspace.getDomain()!=getDomain())
1052     throw DataException("Error - gradient cannot be calculated on different domains.");
1053 jfenwick 1796 DataTypes::ShapeType grad_shape=getDataPointShape();
1054 jgs 94 grad_shape.push_back(functionspace.getDim());
1055     Data out(0.0,grad_shape,functionspace,true);
1056 jfenwick 1872 getDomain()->setToGradient(out,*this);
1057 matt 1353 blocktimer_increment("grad()", blocktimer_start);
1058 jgs 94 return out;
1059     }
1060    
1061     Data
1062     Data::grad() const
1063     {
1064 jfenwick 1803 if (isEmpty())
1065     {
1066     throw DataException("Error - operation not permitted on instances of DataEmpty.");
1067     }
1068 jfenwick 1872 return gradOn(escript::function(*getDomain()));
1069 jgs 94 }
1070    
1071     int
1072     Data::getDataPointSize() const
1073     {
1074 jfenwick 1796 return m_data->getNoValues();
1075 jgs 94 }
1076    
1077 jfenwick 2459
1078 jfenwick 1796 DataTypes::ValueType::size_type
1079 jgs 94 Data::getLength() const
1080     {
1081     return m_data->getLength();
1082     }
1083    
1084 jfenwick 2481
1085 jfenwick 2459 // There is no parallelism here ... elements need to be added in the correct order.
1086     // If we could presize the list and then fill in the elements it might work
1087     // This would need setting elements to be threadsafe.
1088     // Having mulitple C threads calling into one interpreter is aparently a no-no.
1089 jfenwick 2271 const boost::python::object
1090 jfenwick 2459 Data::toListOfTuples(bool scalarastuple)
1091     {
1092     using namespace boost::python;
1093 jfenwick 2461 using boost::python::list;
1094 jfenwick 2459 if (get_MPISize()>1)
1095     {
1096     throw DataException("::toListOfTuples is not available for MPI with more than one process.");
1097     }
1098     unsigned int rank=getDataPointRank();
1099     unsigned int size=getDataPointSize();
1100 jfenwick 2461
1101 jfenwick 2459 int npoints=getNumDataPoints();
1102     expand(); // This will also resolve if required
1103     const DataTypes::ValueType& vec=getReady()->getVectorRO();
1104 jfenwick 2461 boost::python::list temp;
1105     temp.append(object());
1106     boost::python::list res(temp*npoints);// presize the list by the "[None] * npoints" trick
1107 jfenwick 2459 if (rank==0)
1108     {
1109     long count;
1110     if (scalarastuple)
1111     {
1112     for (count=0;count<npoints;++count)
1113     {
1114 jfenwick 2461 res[count]=make_tuple(vec[count]);
1115 jfenwick 2459 }
1116     }
1117     else
1118     {
1119     for (count=0;count<npoints;++count)
1120     {
1121 jfenwick 2461 res[count]=vec[count];
1122 jfenwick 2459 }
1123     }
1124     }
1125     else if (rank==1)
1126     {
1127     size_t count;
1128     size_t offset=0;
1129     for (count=0;count<npoints;++count,offset+=size)
1130     {
1131 jfenwick 2461 res[count]=pointToTuple1(getDataPointShape(), vec, offset);
1132 jfenwick 2459 }
1133     }
1134     else if (rank==2)
1135     {
1136     size_t count;
1137     size_t offset=0;
1138     for (count=0;count<npoints;++count,offset+=size)
1139     {
1140 jfenwick 2461 res[count]=pointToTuple2(getDataPointShape(), vec, offset);
1141 jfenwick 2459 }
1142     }
1143     else if (rank==3)
1144     {
1145     size_t count;
1146     size_t offset=0;
1147     for (count=0;count<npoints;++count,offset+=size)
1148     {
1149 jfenwick 2461 res[count]=pointToTuple3(getDataPointShape(), vec, offset);
1150 jfenwick 2459 }
1151     }
1152     else if (rank==4)
1153     {
1154     size_t count;
1155     size_t offset=0;
1156     for (count=0;count<npoints;++count,offset+=size)
1157     {
1158 jfenwick 2461 res[count]=pointToTuple4(getDataPointShape(), vec, offset);
1159 jfenwick 2459 }
1160     }
1161     else
1162     {
1163     throw DataException("Unknown rank in ::toListOfTuples()");
1164     }
1165     return res;
1166     }
1167    
1168     const boost::python::object
1169 jfenwick 2271 Data::getValueOfDataPointAsTuple(int dataPointNo)
1170     {
1171     forceResolve();
1172     if (getNumDataPointsPerSample()>0) {
1173 gross 921 int sampleNo = dataPointNo/getNumDataPointsPerSample();
1174     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1175     //
1176     // Check a valid sample number has been supplied
1177 trankine 924 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1178 jfenwick 2271 throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.");
1179 gross 921 }
1180 ksteube 1312
1181 gross 921 //
1182     // Check a valid data point number has been supplied
1183 trankine 924 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1184 jfenwick 2271 throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid dataPointNoInSample.");
1185 gross 921 }
1186     // TODO: global error handling
1187 jfenwick 1796 DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1188 jfenwick 2271 return pointToTuple(getDataPointShape(),&(getDataAtOffsetRO(offset)));
1189     }
1190     else
1191     {
1192 jfenwick 2471 // The pre-numpy method would return an empty array of the given shape
1193 jfenwick 2271 // I'm going to throw an exception because if we have zero points per sample we have problems
1194     throw DataException("Error - need at least 1 datapoint per sample.");
1195     }
1196 ksteube 1312
1197 jgs 117 }
1198 jfenwick 1796
1199 jfenwick 2271
1200 gross 921 void
1201 ksteube 1312 Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
1202 jgs 121 {
1203 gross 1034 // this will throw if the value cannot be represented
1204 jfenwick 2271 setValueOfDataPointToArray(dataPointNo,py_object);
1205 gross 1034 }
1206    
1207     void
1208 jfenwick 2271 Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::object& obj)
1209 gross 1034 {
1210 gross 921 if (isProtected()) {
1211     throw DataException("Error - attempt to update protected Data object.");
1212     }
1213 jfenwick 2271 forceResolve();
1214    
1215     WrappedArray w(obj);
1216 gross 921 //
1217     // check rank
1218 jfenwick 2271 if (static_cast<unsigned int>(w.getRank())<getDataPointRank())
1219 jfenwick 2458 throw DataException("Rank of array does not match Data object rank");
1220 bcumming 790
1221 jgs 121 //
1222 jfenwick 2458 // check shape of array
1223 phornby 1953 for (unsigned int i=0; i<getDataPointRank(); i++) {
1224 jfenwick 2271 if (w.getShape()[i]!=getDataPointShape()[i])
1225 jfenwick 2458 throw DataException("Shape of array does not match Data object rank");
1226 jgs 121 }
1227     //
1228 gross 921 // make sure data is expanded:
1229 gross 1859 //
1230 gross 921 if (!isExpanded()) {
1231     expand();
1232 jgs 121 }
1233 gross 921 if (getNumDataPointsPerSample()>0) {
1234     int sampleNo = dataPointNo/getNumDataPointsPerSample();
1235     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1236 jfenwick 2271 m_data->copyToDataPoint(sampleNo, dataPointNoInSample,w);
1237 gross 921 } else {
1238 jfenwick 2271 m_data->copyToDataPoint(-1, 0,w);
1239 gross 921 }
1240     }
1241 jgs 121
1242 gross 922 void
1243     Data::setValueOfDataPoint(int dataPointNo, const double value)
1244     {
1245     if (isProtected()) {
1246     throw DataException("Error - attempt to update protected Data object.");
1247     }
1248     //
1249     // make sure data is expanded:
1250 jfenwick 2271 forceResolve();
1251 gross 922 if (!isExpanded()) {
1252     expand();
1253     }
1254     if (getNumDataPointsPerSample()>0) {
1255     int sampleNo = dataPointNo/getNumDataPointsPerSample();
1256     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1257     m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
1258     } else {
1259     m_data->copyToDataPoint(-1, 0,value);
1260     }
1261     }
1262    
1263 ksteube 1312 const
1264 jfenwick 2271 boost::python::object
1265     Data::getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo)
1266     {
1267     // This could be lazier than it is now
1268     forceResolve();
1269 jgs 121
1270 jfenwick 2271 // copy datapoint into a buffer
1271     // broadcast buffer to all nodes
1272     // convert buffer to tuple
1273     // return tuple
1274 jgs 121
1275 jfenwick 2271 const DataTypes::ShapeType& dataPointShape = getDataPointShape();
1276     size_t length=DataTypes::noValues(dataPointShape);
1277 jgs 121
1278 gross 921 // added for the MPI communication
1279     double *tmpData = new double[length];
1280 bcumming 790
1281 jfenwick 2271 // updated for the MPI case
1282     if( get_MPIRank()==procNo ){
1283     if (getNumDataPointsPerSample()>0) {
1284     int sampleNo = dataPointNo/getNumDataPointsPerSample();
1285     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1286     //
1287     // Check a valid sample number has been supplied
1288     if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1289 jfenwick 2471 throw DataException("Error - Data::getValueOfGlobalDataPointAsTuple: invalid sampleNo.");
1290 bcumming 790 }
1291    
1292 jfenwick 2271 //
1293     // Check a valid data point number has been supplied
1294     if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1295 jfenwick 2471 throw DataException("Error - Data::getValueOfGlobalDataPointAsTuple: invalid dataPointNoInSample.");
1296 bcumming 790 }
1297 jfenwick 2271 // TODO: global error handling
1298     DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1299 bcumming 790
1300 jfenwick 2271 memcpy(tmpData,&(getDataAtOffsetRO(offset)),length*sizeof(double));
1301     }
1302     }
1303     #ifdef PASO_MPI
1304     // broadcast the data to all other processes
1305     MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1306     #endif
1307    
1308     boost::python::tuple t=pointToTuple(dataPointShape,tmpData);
1309     delete [] tmpData;
1310 jgs 121 //
1311     // return the loaded array
1312 jfenwick 2271 return t;
1313    
1314 jgs 121 }
1315    
1316 gross 921
1317 jfenwick 2271 boost::python::object
1318     Data::integrateToTuple_const() const
1319 jfenwick 2005 {
1320     if (isLazy())
1321     {
1322     throw DataException("Error - cannot integrate for constant lazy data.");
1323     }
1324     return integrateWorker();
1325     }
1326 gross 921
1327 jfenwick 2271 boost::python::object
1328     Data::integrateToTuple()
1329 jgs 94 {
1330 jfenwick 2005 if (isLazy())
1331     {
1332     expand();
1333     }
1334     return integrateWorker();
1335 jfenwick 2271
1336 jfenwick 2005 }
1337    
1338 jfenwick 2271 boost::python::object
1339 jfenwick 2005 Data::integrateWorker() const
1340     {
1341 jfenwick 1796 DataTypes::ShapeType shape = getDataPointShape();
1342 ksteube 1312 int dataPointSize = getDataPointSize();
1343 jgs 94
1344     //
1345     // calculate the integral values
1346 ksteube 1312 vector<double> integrals(dataPointSize);
1347     vector<double> integrals_local(dataPointSize);
1348 jfenwick 2089 const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1349     if (dom==0)
1350     {
1351     throw DataException("Can not integrate over non-continuous domains.");
1352     }
1353 ksteube 1312 #ifdef PASO_MPI
1354 jfenwick 2089 dom->setToIntegrals(integrals_local,*this);
1355 ksteube 1312 // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1356     double *tmp = new double[dataPointSize];
1357     double *tmp_local = new double[dataPointSize];
1358     for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1359     MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1360     for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1361 jfenwick 2271 tuple result=pointToTuple(shape,tmp);
1362 ksteube 1312 delete[] tmp;
1363     delete[] tmp_local;
1364     #else
1365 jfenwick 2089 dom->setToIntegrals(integrals,*this);
1366 jfenwick 2271 /* double *tmp = new double[dataPointSize];
1367     for (int i=0; i<dataPointSize; i++) { tmp[i]=integrals[i]; }*/
1368     tuple result=pointToTuple(shape,integrals);
1369     // delete tmp;
1370 ksteube 1312 #endif
1371 jgs 94
1372    
1373 jfenwick 2271 return result;
1374 jgs 94 }
1375    
1376     Data
1377     Data::sin() const
1378     {
1379 jfenwick 2146 MAKELAZYOP(SIN)
1380 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1381 jgs 94 }
1382    
1383     Data
1384     Data::cos() const
1385     {
1386 jfenwick 2146 MAKELAZYOP(COS)
1387 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1388 jgs 94 }
1389    
1390     Data
1391     Data::tan() const
1392     {
1393 jfenwick 2146 MAKELAZYOP(TAN)
1394 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1395 jgs 94 }
1396    
1397     Data
1398 jgs 150 Data::asin() const
1399     {
1400 jfenwick 2146 MAKELAZYOP(ASIN)
1401 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1402 jgs 150 }
1403    
1404     Data
1405     Data::acos() const
1406     {
1407 jfenwick 2146 MAKELAZYOP(ACOS)
1408 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1409 jgs 150 }
1410    
1411 phornby 1026
1412 jgs 150 Data
1413     Data::atan() const
1414     {
1415 jfenwick 2146 MAKELAZYOP(ATAN)
1416 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1417 jgs 150 }
1418    
1419     Data
1420     Data::sinh() const
1421     {
1422 jfenwick 2146 MAKELAZYOP(SINH)
1423 jfenwick 2005 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1424 jgs 150 }
1425    
1426     Data
1427     Data::cosh() const
1428     {
1429 jfenwick 2146 MAKELAZYOP(COSH)
1430 jfenwick 2005 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1431 jgs 150 }
1432    
1433     Data
1434     Data::tanh() const
1435     {
1436 jfenwick 2146 MAKELAZYOP(TANH)
1437 jfenwick 2005 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1438 jgs 150 }
1439    
1440 phornby 1026
1441 jgs 150 Data
1442 phornby 1026 Data::erf() const
1443     {
1444 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1445 gross 1028 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1446     #else
1447 jfenwick 2146 MAKELAZYOP(ERF)
1448 matt 1334 return C_TensorUnaryOperation(*this, ::erf);
1449 phornby 1026 #endif
1450     }
1451    
1452     Data
1453 jgs 150 Data::asinh() const
1454     {
1455 jfenwick 2146 MAKELAZYOP(ASINH)
1456 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1457 matt 1334 return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1458 phornby 1032 #else
1459 matt 1334 return C_TensorUnaryOperation(*this, ::asinh);
1460 phornby 1032 #endif
1461 jgs 150 }
1462    
1463     Data
1464     Data::acosh() const
1465     {
1466 jfenwick 2146 MAKELAZYOP(ACOSH)
1467 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1468 matt 1334 return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1469 phornby 1032 #else
1470 matt 1334 return C_TensorUnaryOperation(*this, ::acosh);
1471 phornby 1032 #endif
1472 jgs 150 }
1473    
1474     Data
1475     Data::atanh() const
1476     {
1477 jfenwick 2146 MAKELAZYOP(ATANH)
1478 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1479 matt 1334 return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1480 phornby 1032 #else
1481 matt 1334 return C_TensorUnaryOperation(*this, ::atanh);
1482 phornby 1032 #endif
1483 jgs 150 }
1484    
1485     Data
1486 gross 286 Data::log10() const
1487 jfenwick 2146 {
1488     MAKELAZYOP(LOG10)
1489 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1490 jgs 94 }
1491    
1492     Data
1493 gross 286 Data::log() const
1494 jgs 94 {
1495 jfenwick 2146 MAKELAZYOP(LOG)
1496 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1497 jgs 94 }
1498    
1499 jgs 106 Data
1500     Data::sign() const
1501 jgs 94 {
1502 jfenwick 2146 MAKELAZYOP(SIGN)
1503 matt 1334 return C_TensorUnaryOperation(*this, escript::fsign);
1504 jgs 94 }
1505    
1506 jgs 106 Data
1507     Data::abs() const
1508 jgs 94 {
1509 jfenwick 2146 MAKELAZYOP(ABS)
1510 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1511 jgs 94 }
1512    
1513 jgs 106 Data
1514     Data::neg() const
1515 jgs 94 {
1516 jfenwick 2146 MAKELAZYOP(NEG)
1517 matt 1334 return C_TensorUnaryOperation(*this, negate<double>());
1518 jgs 94 }
1519    
1520 jgs 102 Data
1521 jgs 106 Data::pos() const
1522 jgs 94 {
1523 jfenwick 2005 // not doing lazy check here is deliberate.
1524     // since a deep copy of lazy data should be cheap, I'll just let it happen now
1525 jgs 148 Data result;
1526     // perform a deep copy
1527     result.copy(*this);
1528     return result;
1529 jgs 102 }
1530    
1531     Data
1532 jgs 106 Data::exp() const
1533 jfenwick 2146 {
1534     MAKELAZYOP(EXP)
1535 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1536 jgs 102 }
1537    
1538     Data
1539 jgs 106 Data::sqrt() const
1540 jgs 102 {
1541 jfenwick 2146 MAKELAZYOP(SQRT)
1542 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1543 jgs 102 }
1544    
1545 jgs 106 double
1546 jfenwick 2005 Data::Lsup_const() const
1547 jgs 102 {
1548 jfenwick 2005 if (isLazy())
1549     {
1550     throw DataException("Error - cannot compute Lsup for constant lazy data.");
1551     }
1552     return LsupWorker();
1553     }
1554    
1555     double
1556     Data::Lsup()
1557     {
1558     if (isLazy())
1559     {
1560 jfenwick 2085 resolve();
1561 jfenwick 2005 }
1562     return LsupWorker();
1563     }
1564    
1565     double
1566     Data::sup_const() const
1567     {
1568     if (isLazy())
1569     {
1570     throw DataException("Error - cannot compute sup for constant lazy data.");
1571     }
1572     return supWorker();
1573     }
1574    
1575     double
1576     Data::sup()
1577     {
1578     if (isLazy())
1579     {
1580 jfenwick 2085 resolve();
1581 jfenwick 2005 }
1582     return supWorker();
1583     }
1584    
1585     double
1586     Data::inf_const() const
1587     {
1588     if (isLazy())
1589     {
1590     throw DataException("Error - cannot compute inf for constant lazy data.");
1591     }
1592     return infWorker();
1593     }
1594    
1595     double
1596     Data::inf()
1597     {
1598     if (isLazy())
1599     {
1600 jfenwick 2085 resolve();
1601 jfenwick 2005 }
1602     return infWorker();
1603     }
1604    
1605     double
1606     Data::LsupWorker() const
1607     {
1608 phornby 1455 double localValue;
1609 jgs 106 //
1610     // set the initial absolute maximum value to zero
1611 bcumming 751
1612 jgs 147 AbsMax abs_max_func;
1613 bcumming 751 localValue = algorithm(abs_max_func,0);
1614     #ifdef PASO_MPI
1615 phornby 1455 double globalValue;
1616 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1617     return globalValue;
1618     #else
1619     return localValue;
1620     #endif
1621 jgs 117 }
1622    
1623     double
1624 jfenwick 2005 Data::supWorker() const
1625 jgs 102 {
1626 phornby 1455 double localValue;
1627 jgs 106 //
1628     // set the initial maximum value to min possible double
1629 jgs 147 FMax fmax_func;
1630 bcumming 751 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1631     #ifdef PASO_MPI
1632 phornby 1455 double globalValue;
1633 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1634     return globalValue;
1635     #else
1636     return localValue;
1637     #endif
1638 jgs 102 }
1639    
1640 jgs 106 double
1641 jfenwick 2005 Data::infWorker() const
1642 jgs 102 {
1643 phornby 1455 double localValue;
1644 jgs 106 //
1645     // set the initial minimum value to max possible double
1646 jgs 147 FMin fmin_func;
1647 bcumming 751 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1648     #ifdef PASO_MPI
1649 phornby 1455 double globalValue;
1650 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1651     return globalValue;
1652     #else
1653     return localValue;
1654     #endif
1655 jgs 102 }
1656    
1657 bcumming 751 /* TODO */
1658     /* global reduction */
1659 jgs 102 Data
1660 jgs 106 Data::maxval() const
1661 jgs 102 {
1662 jfenwick 2037 if (isLazy())
1663     {
1664     Data temp(*this); // to get around the fact that you can't resolve a const Data
1665     temp.resolve();
1666     return temp.maxval();
1667     }
1668 jgs 113 //
1669     // set the initial maximum value to min possible double
1670 jgs 147 FMax fmax_func;
1671     return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1672 jgs 102 }
1673    
1674     Data
1675 jgs 106 Data::minval() const
1676 jgs 102 {
1677 jfenwick 2037 if (isLazy())
1678     {
1679     Data temp(*this); // to get around the fact that you can't resolve a const Data
1680     temp.resolve();
1681     return temp.minval();
1682     }
1683 jgs 113 //
1684     // set the initial minimum value to max possible double
1685 jgs 147 FMin fmin_func;
1686     return dp_algorithm(fmin_func,numeric_limits<double>::max());
1687 jgs 102 }
1688    
1689 jgs 123 Data
1690 gross 804 Data::swapaxes(const int axis0, const int axis1) const
1691 jgs 123 {
1692 gross 804 int axis0_tmp,axis1_tmp;
1693 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1694     DataTypes::ShapeType ev_shape;
1695 gross 800 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1696     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1697     int rank=getDataPointRank();
1698 gross 804 if (rank<2) {
1699     throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1700 gross 800 }
1701 gross 804 if (axis0<0 || axis0>rank-1) {
1702     throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1703     }
1704     if (axis1<0 || axis1>rank-1) {
1705     throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1706     }
1707     if (axis0 == axis1) {
1708     throw DataException("Error - Data::swapaxes: axis indices must be different.");
1709     }
1710 jfenwick 2496 MAKELAZYOP2(SWAP,axis0,axis1)
1711     if (axis0 > axis1)
1712     {
1713     axis0_tmp=axis1;
1714     axis1_tmp=axis0;
1715 gross 804 }
1716 jfenwick 2496 else
1717     {
1718     axis0_tmp=axis0;
1719     axis1_tmp=axis1;
1720 gross 800 }
1721 jfenwick 2496 for (int i=0; i<rank; i++)
1722     {
1723     if (i == axis0_tmp)
1724     {
1725     ev_shape.push_back(s[axis1_tmp]);
1726     }
1727     else if (i == axis1_tmp)
1728     {
1729     ev_shape.push_back(s[axis0_tmp]);
1730     }
1731     else
1732     {
1733     ev_shape.push_back(s[i]);
1734     }
1735     }
1736 gross 800 Data ev(0.,ev_shape,getFunctionSpace());
1737     ev.typeMatchRight(*this);
1738 gross 804 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1739 gross 800 return ev;
1740 jgs 123 }
1741    
1742     Data
1743 ksteube 775 Data::symmetric() const
1744 jgs 123 {
1745 ksteube 775 // check input
1746 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1747 ksteube 775 if (getDataPointRank()==2) {
1748 ksteube 1312 if(s[0] != s[1])
1749 ksteube 775 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1750     }
1751     else if (getDataPointRank()==4) {
1752     if(!(s[0] == s[2] && s[1] == s[3]))
1753     throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1754     }
1755     else {
1756     throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1757     }
1758 jfenwick 2146 MAKELAZYOP(SYM)
1759 ksteube 775 Data ev(0.,getDataPointShape(),getFunctionSpace());
1760     ev.typeMatchRight(*this);
1761     m_data->symmetric(ev.m_data.get());
1762     return ev;
1763     }
1764    
1765     Data
1766     Data::nonsymmetric() const
1767     {
1768 jfenwick 2146 MAKELAZYOP(NSYM)
1769 ksteube 775 // check input
1770 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1771 ksteube 775 if (getDataPointRank()==2) {
1772 ksteube 1312 if(s[0] != s[1])
1773 ksteube 775 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1774 jfenwick 1796 DataTypes::ShapeType ev_shape;
1775 ksteube 775 ev_shape.push_back(s[0]);
1776     ev_shape.push_back(s[1]);
1777     Data ev(0.,ev_shape,getFunctionSpace());
1778     ev.typeMatchRight(*this);
1779     m_data->nonsymmetric(ev.m_data.get());
1780     return ev;
1781     }
1782     else if (getDataPointRank()==4) {
1783     if(!(s[0] == s[2] && s[1] == s[3]))
1784     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1785 jfenwick 1796 DataTypes::ShapeType ev_shape;
1786 ksteube 775 ev_shape.push_back(s[0]);
1787     ev_shape.push_back(s[1]);
1788     ev_shape.push_back(s[2]);
1789     ev_shape.push_back(s[3]);
1790     Data ev(0.,ev_shape,getFunctionSpace());
1791     ev.typeMatchRight(*this);
1792     m_data->nonsymmetric(ev.m_data.get());
1793     return ev;
1794     }
1795     else {
1796     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1797     }
1798     }
1799    
1800     Data
1801 gross 800 Data::trace(int axis_offset) const
1802 jfenwick 2084 {
1803 jfenwick 2146 MAKELAZYOPOFF(TRACE,axis_offset)
1804 jfenwick 2200 if ((axis_offset<0) || (axis_offset>getDataPointRank()))
1805     {
1806     throw DataException("Error - Data::trace, axis_offset must be between 0 and rank-2 inclusive.");
1807     }
1808 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1809 ksteube 775 if (getDataPointRank()==2) {
1810 jfenwick 1796 DataTypes::ShapeType ev_shape;
1811 ksteube 775 Data ev(0.,ev_shape,getFunctionSpace());
1812     ev.typeMatchRight(*this);
1813 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1814 ksteube 775 return ev;
1815     }
1816     if (getDataPointRank()==3) {
1817 jfenwick 1796 DataTypes::ShapeType ev_shape;
1818 ksteube 775 if (axis_offset==0) {
1819     int s2=s[2];
1820     ev_shape.push_back(s2);
1821     }
1822     else if (axis_offset==1) {
1823     int s0=s[0];
1824     ev_shape.push_back(s0);
1825     }
1826     Data ev(0.,ev_shape,getFunctionSpace());
1827     ev.typeMatchRight(*this);
1828 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1829 ksteube 775 return ev;
1830     }
1831     if (getDataPointRank()==4) {
1832 jfenwick 1796 DataTypes::ShapeType ev_shape;
1833 ksteube 775 if (axis_offset==0) {
1834     ev_shape.push_back(s[2]);
1835     ev_shape.push_back(s[3]);
1836     }
1837     else if (axis_offset==1) {
1838     ev_shape.push_back(s[0]);
1839     ev_shape.push_back(s[3]);
1840     }
1841     else if (axis_offset==2) {
1842     ev_shape.push_back(s[0]);
1843     ev_shape.push_back(s[1]);
1844     }
1845     Data ev(0.,ev_shape,getFunctionSpace());
1846     ev.typeMatchRight(*this);
1847 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1848 ksteube 775 return ev;
1849     }
1850     else {
1851 gross 800 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1852 ksteube 775 }
1853     }
1854    
1855     Data
1856     Data::transpose(int axis_offset) const
1857 jfenwick 2005 {
1858 jfenwick 2146 MAKELAZYOPOFF(TRANS,axis_offset)
1859 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1860     DataTypes::ShapeType ev_shape;
1861 ksteube 775 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1862     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1863     int rank=getDataPointRank();
1864     if (axis_offset<0 || axis_offset>rank) {
1865     throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1866     }
1867     for (int i=0; i<rank; i++) {
1868     int index = (axis_offset+i)%rank;
1869     ev_shape.push_back(s[index]); // Append to new shape
1870     }
1871     Data ev(0.,ev_shape,getFunctionSpace());
1872     ev.typeMatchRight(*this);
1873     m_data->transpose(ev.m_data.get(), axis_offset);
1874     return ev;
1875 jgs 123 }
1876    
1877 gross 576 Data
1878     Data::eigenvalues() const
1879     {
1880 jfenwick 2005 if (isLazy())
1881     {
1882     Data temp(*this); // to get around the fact that you can't resolve a const Data
1883     temp.resolve();
1884     return temp.eigenvalues();
1885     }
1886 gross 576 // check input
1887 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1888 ksteube 1312 if (getDataPointRank()!=2)
1889 gross 576 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1890 ksteube 1312 if(s[0] != s[1])
1891 gross 576 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1892     // create return
1893 jfenwick 1796 DataTypes::ShapeType ev_shape(1,s[0]);
1894 gross 576 Data ev(0.,ev_shape,getFunctionSpace());
1895     ev.typeMatchRight(*this);
1896     m_data->eigenvalues(ev.m_data.get());
1897     return ev;
1898     }
1899    
1900 jgs 121 const boost::python::tuple
1901 gross 576 Data::eigenvalues_and_eigenvectors(const double tol) const
1902     {
1903 jfenwick 2005 if (isLazy())
1904     {
1905     Data temp(*this); // to get around the fact that you can't resolve a const Data
1906     temp.resolve();
1907     return temp.eigenvalues_and_eigenvectors(tol);
1908     }
1909 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1910 ksteube 1312 if (getDataPointRank()!=2)
1911 gross 576 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1912 ksteube 1312 if(s[0] != s[1])
1913 gross 576 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1914     // create return
1915 jfenwick 1796 DataTypes::ShapeType ev_shape(1,s[0]);
1916 gross 576 Data ev(0.,ev_shape,getFunctionSpace());
1917     ev.typeMatchRight(*this);
1918 jfenwick 1796 DataTypes::ShapeType V_shape(2,s[0]);
1919 gross 576 Data V(0.,V_shape,getFunctionSpace());
1920     V.typeMatchRight(*this);
1921     m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1922     return make_tuple(boost::python::object(ev),boost::python::object(V));
1923     }
1924    
1925     const boost::python::tuple
1926 gross 921 Data::minGlobalDataPoint() const
1927 jgs 121 {
1928 gross 921 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1929 jgs 148 // abort (for unknown reasons) if there are openmp directives with it in the
1930     // surrounding function
1931    
1932     int DataPointNo;
1933 gross 921 int ProcNo;
1934     calc_minGlobalDataPoint(ProcNo,DataPointNo);
1935     return make_tuple(ProcNo,DataPointNo);
1936 jgs 148 }
1937    
1938     void
1939 gross 921 Data::calc_minGlobalDataPoint(int& ProcNo,
1940     int& DataPointNo) const
1941 jgs 148 {
1942 jfenwick 2005 if (isLazy())
1943     {
1944     Data temp(*this); // to get around the fact that you can't resolve a const Data
1945     temp.resolve();
1946     return temp.calc_minGlobalDataPoint(ProcNo,DataPointNo);
1947     }
1948 jgs 148 int i,j;
1949     int lowi=0,lowj=0;
1950     double min=numeric_limits<double>::max();
1951    
1952 jgs 121 Data temp=minval();
1953    
1954     int numSamples=temp.getNumSamples();
1955     int numDPPSample=temp.getNumDataPointsPerSample();
1956    
1957 jgs 148 double next,local_min;
1958 jfenwick 1946 int local_lowi=0,local_lowj=0;
1959 jgs 121
1960 caltinay 2081 #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
1961 jgs 148 {
1962     local_min=min;
1963     #pragma omp for private(i,j) schedule(static)
1964     for (i=0; i<numSamples; i++) {
1965     for (j=0; j<numDPPSample; j++) {
1966 jfenwick 2271 next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
1967 jgs 148 if (next<local_min) {
1968     local_min=next;
1969     local_lowi=i;
1970     local_lowj=j;
1971     }
1972 jgs 121 }
1973     }
1974 jgs 148 #pragma omp critical
1975 jfenwick 2476 if (local_min<min) { // If we found a smaller value than our sentinel
1976 jgs 148 min=local_min;
1977     lowi=local_lowi;
1978     lowj=local_lowj;
1979     }
1980 jgs 121 }
1981    
1982 bcumming 782 #ifdef PASO_MPI
1983 jfenwick 2476 // determine the processor on which the minimum occurs
1984     next = temp.getDataPointRO(lowi,lowj);
1985     int lowProc = 0;
1986     double *globalMins = new double[get_MPISize()+1];
1987     int error;
1988     error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1989 ksteube 1312
1990 jfenwick 2476 if( get_MPIRank()==0 ){
1991     next = globalMins[lowProc];
1992     for( i=1; i<get_MPISize(); i++ )
1993     if( next>globalMins[i] ){
1994     lowProc = i;
1995     next = globalMins[i];
1996     }
1997     }
1998 gross 2492 MPI_Bcast( &lowProc, 1, MPI_INT, 0, get_MPIComm() );
1999 bcumming 782
2000 jfenwick 2476 delete [] globalMins;
2001     ProcNo = lowProc;
2002 bcumming 790 #else
2003 jfenwick 2476 ProcNo = 0;
2004 bcumming 782 #endif
2005 gross 921 DataPointNo = lowj + lowi * numDPPSample;
2006 jgs 121 }
2007    
2008 jfenwick 2476
2009     const boost::python::tuple
2010     Data::maxGlobalDataPoint() const
2011     {
2012     int DataPointNo;
2013     int ProcNo;
2014     calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2015     return make_tuple(ProcNo,DataPointNo);
2016     }
2017    
2018 jgs 104 void
2019 jfenwick 2476 Data::calc_maxGlobalDataPoint(int& ProcNo,
2020     int& DataPointNo) const
2021     {
2022     if (isLazy())
2023     {
2024     Data temp(*this); // to get around the fact that you can't resolve a const Data
2025     temp.resolve();
2026     return temp.calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2027     }
2028     int i,j;
2029     int highi=0,highj=0;
2030     //-------------
2031     double max=numeric_limits<double>::min();
2032    
2033     Data temp=maxval();
2034    
2035     int numSamples=temp.getNumSamples();
2036     int numDPPSample=temp.getNumDataPointsPerSample();
2037    
2038     double next,local_max;
2039     int local_highi=0,local_highj=0;
2040    
2041     #pragma omp parallel firstprivate(local_highi,local_highj) private(next,local_max)
2042     {
2043     local_max=max;
2044     #pragma omp for private(i,j) schedule(static)
2045     for (i=0; i<numSamples; i++) {
2046     for (j=0; j<numDPPSample; j++) {
2047     next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2048     if (next>local_max) {
2049     local_max=next;
2050     local_highi=i;
2051     local_highj=j;
2052     }
2053     }
2054     }
2055     #pragma omp critical
2056     if (local_max>max) { // If we found a larger value than our sentinel
2057     max=local_max;
2058     highi=local_highi;
2059     highj=local_highj;
2060     }
2061     }
2062    
2063     #ifdef PASO_MPI
2064     // determine the processor on which the maximum occurs
2065     next = temp.getDataPointRO(highi,highj);
2066     int highProc = 0;
2067     double *globalMaxs = new double[get_MPISize()+1];
2068     int error;
2069     error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMaxs, 1, MPI_DOUBLE, 0, get_MPIComm() );
2070    
2071     if( get_MPIRank()==0 ){
2072     next = globalMaxs[highProc];
2073     for( i=1; i<get_MPISize(); i++ )
2074     if( next>globalMaxs[i] ){
2075     highProc = i;
2076     next = globalMaxs[i];
2077     }
2078     }
2079 gross 2492 MPI_Bcast( &highProc, 1, MPI_INT, 0, get_MPIComm() );
2080 jfenwick 2476 delete [] globalMaxs;
2081     ProcNo = highProc;
2082     #else
2083     ProcNo = 0;
2084     #endif
2085     DataPointNo = highj + highi * numDPPSample;
2086     }
2087    
2088     void
2089 jgs 104 Data::saveDX(std::string fileName) const
2090     {
2091 jfenwick 1803 if (isEmpty())
2092     {
2093     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2094     }
2095 jfenwick 2005 if (isLazy())
2096     {
2097     Data temp(*this); // to get around the fact that you can't resolve a const Data
2098     temp.resolve();
2099     temp.saveDX(fileName);
2100     return;
2101     }
2102 jgs 153 boost::python::dict args;
2103     args["data"]=boost::python::object(this);
2104 jfenwick 1872 getDomain()->saveDX(fileName,args);
2105 jgs 104 return;
2106     }
2107    
2108 jgs 110 void
2109     Data::saveVTK(std::string fileName) const
2110     {
2111 jfenwick 1803 if (isEmpty())
2112     {
2113     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2114     }
2115 jfenwick 2005 if (isLazy())
2116     {
2117     Data temp(*this); // to get around the fact that you can't resolve a const Data
2118     temp.resolve();
2119     temp.saveVTK(fileName);
2120     return;
2121     }
2122 jgs 153 boost::python::dict args;
2123     args["data"]=boost::python::object(this);
2124 gross 2421 getDomain()->saveVTK(fileName,args,"","");
2125 jgs 110 return;
2126     }
2127    
2128 jfenwick 2635
2129    
2130 jgs 102 Data&
2131     Data::operator+=(const Data& right)
2132     {
2133 gross 783 if (isProtected()) {
2134     throw DataException("Error - attempt to update protected Data object.");
2135     }
2136 jfenwick 2146 MAKELAZYBINSELF(right,ADD) // for lazy + is equivalent to +=
2137 jfenwick 2271 exclusiveWrite(); // Since Lazy data does not modify its leaves we only need to worry here
2138 jfenwick 2146 binaryOp(right,plus<double>());
2139     return (*this);
2140 jgs 94 }
2141    
2142 jgs 102 Data&
2143     Data::operator+=(const boost::python::object& right)
2144 jgs 94 {
2145 jfenwick 2005 if (isProtected()) {
2146     throw DataException("Error - attempt to update protected Data object.");
2147     }
2148 gross 854 Data tmp(right,getFunctionSpace(),false);
2149 jfenwick 2271 (*this)+=tmp;
2150     return *this;
2151 jgs 94 }
2152 jfenwick 2005
2153     // Hmmm, operator= makes a deep copy but the copy constructor does not?
2154 ksteube 1312 Data&
2155     Data::operator=(const Data& other)
2156     {
2157 jfenwick 2271 #if defined ASSIGNMENT_MEANS_DEEPCOPY
2158     // This should not be used.
2159     // Just leaving this here until I have completed transition
2160 ksteube 1312 copy(other);
2161 jfenwick 2271 #else
2162     m_protected=false; // since any changes should be caught by exclusiveWrite();
2163     // m_data=other.m_data;
2164     set_m_data(other.m_data);
2165     #endif
2166 ksteube 1312 return (*this);
2167     }
2168 jgs 94
2169 jgs 102 Data&
2170     Data::operator-=(const Data& right)
2171 jgs 94 {
2172 gross 783 if (isProtected()) {
2173     throw DataException("Error - attempt to update protected Data object.");
2174     }
2175 jfenwick 2146 MAKELAZYBINSELF(right,SUB)
2176 jfenwick 2271 exclusiveWrite();
2177 jfenwick 2146 binaryOp(right,minus<double>());
2178     return (*this);
2179 jgs 94 }
2180    
2181 jgs 102 Data&
2182     Data::operator-=(const boost::python::object& right)
2183 jgs 94 {
2184 jfenwick 2005 if (isProtected()) {
2185     throw DataException("Error - attempt to update protected Data object.");
2186     }
2187 gross 854 Data tmp(right,getFunctionSpace(),false);
2188 jfenwick 2271 (*this)-=tmp;
2189 jfenwick 2146 return (*this);
2190 jgs 94 }
2191    
2192 jgs 102 Data&
2193     Data::operator*=(const Data& right)
2194 jgs 94 {
2195 gross 783 if (isProtected()) {
2196     throw DataException("Error - attempt to update protected Data object.");
2197     }
2198 jfenwick 2146 MAKELAZYBINSELF(right,MUL)
2199 jfenwick 2271 exclusiveWrite();
2200 jfenwick 2146 binaryOp(right,multiplies<double>());
2201     return (*this);
2202 jgs 94 }
2203    
2204 jgs 102 Data&
2205     Data::operator*=(const boost::python::object& right)
2206 jfenwick 2005 {
2207     if (isProtected()) {
2208     throw DataException("Error - attempt to update protected Data object.");
2209     }
2210 gross 854 Data tmp(right,getFunctionSpace(),false);
2211 jfenwick 2271 (*this)*=tmp;
2212 jfenwick 2146 return (*this);
2213 jgs 94 }
2214    
2215 jgs 102 Data&
2216     Data::operator/=(const Data& right)
2217 jgs 94 {
2218 gross 783 if (isProtected()) {
2219     throw DataException("Error - attempt to update protected Data object.");
2220     }
2221 jfenwick 2146 MAKELAZYBINSELF(right,DIV)
2222 jfenwick 2271 exclusiveWrite();
2223 jfenwick 2146 binaryOp(right,divides<double>());
2224     return (*this);
2225 jgs 94 }
2226    
2227 jgs 102 Data&
2228     Data::operator/=(const boost::python::object& right)
2229 jgs 94 {
2230 jfenwick 2005 if (isProtected()) {
2231     throw DataException("Error - attempt to update protected Data object.");
2232     }
2233 gross 854 Data tmp(right,getFunctionSpace(),false);
2234 jfenwick 2271 (*this)/=tmp;
2235 jfenwick 2146 return (*this);
2236 jgs 94 }
2237    
2238 jgs 102 Data
2239 gross 699 Data::rpowO(const boost::python::object& left) const
2240     {
2241     Data left_d(left,*this);
2242     return left_d.powD(*this);
2243     }
2244    
2245     Data
2246 jgs 102 Data::powO(const boost::python::object& right) const
2247 jgs 94 {
2248 gross 854 Data tmp(right,getFunctionSpace(),false);
2249     return powD(tmp);
2250 jgs 94 }
2251    
2252 jgs 102 Data
2253     Data::powD(const Data& right) const
2254 jgs 94 {
2255 jfenwick 2146 MAKELAZYBIN(right,POW)
2256 matt 1350 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2257 jgs 94 }
2258    
2259     //
2260 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2261 jgs 102 Data
2262     escript::operator+(const Data& left, const Data& right)
2263 jgs 94 {
2264 jfenwick 2146 MAKELAZYBIN2(left,right,ADD)
2265 matt 1327 return C_TensorBinaryOperation(left, right, plus<double>());
2266 jgs 94 }
2267    
2268     //
2269 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2270 jgs 102 Data
2271     escript::operator-(const Data& left, const Data& right)
2272 jgs 94 {
2273 jfenwick 2146 MAKELAZYBIN2(left,right,SUB)
2274 matt 1327 return C_TensorBinaryOperation(left, right, minus<double>());
2275 jgs 94 }
2276    
2277     //
2278 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2279 jgs 102 Data
2280     escript::operator*(const Data& left, const Data& right)
2281 jgs 94 {
2282 jfenwick 2146 MAKELAZYBIN2(left,right,MUL)
2283 matt 1327 return C_TensorBinaryOperation(left, right, multiplies<double>());
2284 jgs 94 }
2285    
2286     //
2287 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2288 jgs 102 Data
2289     escript::operator/(const Data& left, const Data& right)
2290 jgs 94 {
2291 jfenwick 2146 MAKELAZYBIN2(left,right,DIV)
2292 matt 1327 return C_TensorBinaryOperation(left, right, divides<double>());
2293 jgs 94 }
2294    
2295     //
2296 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2297 jgs 102 Data
2298     escript::operator+(const Data& left, const boost::python::object& right)
2299 jgs 94 {
2300 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2301     MAKELAZYBIN2(left,tmp,ADD)
2302     return left+tmp;
2303 jgs 94 }
2304    
2305     //
2306 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2307 jgs 102 Data
2308     escript::operator-(const Data& left, const boost::python::object& right)
2309 jgs 94 {
2310 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2311     MAKELAZYBIN2(left,tmp,SUB)
2312     return left-tmp;
2313 jgs 94 }
2314    
2315     //
2316 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2317 jgs 102 Data
2318     escript::operator*(const Data& left, const boost::python::object& right)
2319 jgs 94 {
2320 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2321     MAKELAZYBIN2(left,tmp,MUL)
2322     return left*tmp;
2323 jgs 94 }
2324    
2325     //
2326 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2327 jgs 102 Data
2328     escript::operator/(const Data& left, const boost::python::object& right)
2329 jgs 94 {
2330 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2331     MAKELAZYBIN2(left,tmp,DIV)
2332     return left/tmp;
2333 jgs 94 }
2334    
2335     //
2336 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2337 jgs 102 Data
2338     escript::operator+(const boost::python::object& left, const Data& right)
2339 jgs 94 {
2340 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2341     MAKELAZYBIN2(tmp,right,ADD)
2342     return tmp+right;
2343 jgs 94 }
2344    
2345     //
2346 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2347 jgs 102 Data
2348     escript::operator-(const boost::python::object& left, const Data& right)
2349 jgs 94 {
2350 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2351     MAKELAZYBIN2(tmp,right,SUB)
2352     return tmp-right;
2353 jgs 94 }
2354    
2355     //
2356 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2357 jgs 102 Data
2358     escript::operator*(const boost::python::object& left, const Data& right)
2359 jgs 94 {
2360 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2361     MAKELAZYBIN2(tmp,right,MUL)
2362     return tmp*right;
2363 jgs 94 }
2364    
2365     //
2366 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2367 jgs 102 Data
2368     escript::operator/(const boost::python::object& left, const Data& right)
2369 jgs 94 {
2370 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2371     MAKELAZYBIN2(tmp,right,DIV)
2372     return tmp/right;
2373 jgs 94 }
2374    
2375 jgs 102
2376 bcumming 751 /* TODO */
2377     /* global reduction */
2378 jgs 102 Data
2379 ksteube 1312 Data::getItem(const boost::python::object& key) const
2380 jgs 94 {
2381    
2382 jfenwick 1796 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2383 jgs 94
2384 jfenwick 1796 if (slice_region.size()!=getDataPointRank()) {
2385 jgs 102 throw DataException("Error - slice size does not match Data rank.");
2386 jgs 94 }
2387    
2388 jgs 102 return getSlice(slice_region);
2389 jgs 94 }
2390    
2391 bcumming 751 /* TODO */
2392     /* global reduction */
2393 jgs 94 Data
2394 jfenwick 1796 Data::getSlice(const DataTypes::RegionType& region) const
2395 jgs 94 {
2396 jgs 102 return Data(*this,region);
2397 jgs 94 }
2398    
2399 bcumming 751 /* TODO */
2400     /* global reduction */
2401 jgs 94 void
2402 jgs 102 Data::setItemO(const boost::python::object& key,
2403     const boost::python::object& value)
2404 jgs 94 {
2405 jgs 102 Data tempData(value,getFunctionSpace());
2406     setItemD(key,tempData);
2407     }
2408    
2409     void
2410     Data::setItemD(const boost::python::object& key,
2411     const Data& value)
2412     {
2413 jfenwick 1796 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2414     if (slice_region.size()!=getDataPointRank()) {
2415 jgs 94 throw DataException("Error - slice size does not match Data rank.");
2416     }
2417 jfenwick 2271 exclusiveWrite();
2418 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
2419     setSlice(Data(value,getFunctionSpace()),slice_region);
2420     } else {
2421     setSlice(value,slice_region);
2422     }
2423 jgs 94 }
2424    
2425     void
2426 jgs 102 Data::setSlice(const Data& value,
2427 jfenwick 1796 const DataTypes::RegionType& region)
2428 jgs 94 {
2429 gross 783 if (isProtected()) {
2430     throw DataException("Error - attempt to update protected Data object.");
2431     }
2432 jfenwick 2271 forceResolve();
2433     exclusiveWrite(); // In case someone finds a way to call this without going through setItemD
2434 jgs 102 Data tempValue(value);
2435     typeMatchLeft(tempValue);
2436     typeMatchRight(tempValue);
2437 jfenwick 2005 getReady()->setSlice(tempValue.m_data.get(),region);
2438 jgs 102 }
2439  <