/[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 1910 - (hide annotations)
Thu Oct 23 03:05:28 2008 UTC (11 years, 9 months ago) by jfenwick
File size: 82102 byte(s)
Branch commit.
Support for ** added.

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