/[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 1943 - (hide annotations)
Wed Oct 29 04:05:14 2008 UTC (10 years, 7 months ago) by jfenwick
File size: 83824 byte(s)
Branch commit
Added interpolation.

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