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

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

Parent Directory Parent Directory | Revision Log Revision Log


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