/[escript]/branches/schroedinger/escript/src/Data.cpp
ViewVC logotype

Annotation of /branches/schroedinger/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1888 - (hide annotations)
Wed Oct 15 04:00:21 2008 UTC (13 years, 1 month ago) by jfenwick
File size: 77927 byte(s)
Branch commit.
Added more binary ops.

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