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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2086 - (hide annotations)
Mon Nov 24 02:38:50 2008 UTC (11 years, 4 months ago) by jfenwick
File size: 85194 byte(s)
Added checks in C_GeneralTensorProduct (Data:: and Delayed forms) as 
well as the DataAbstract Constructor to prevent Objects with Rank>4 
being created.

Moved the relevant #define into systemdep.

Removed some comments.

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 1796 throw FunctionSpaceException("Call to probeInterpolation returned false for 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     temp_data.copyFromNumArray(asNumArray);
215    
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     Data*
411     Data::copySelf()
412     {
413     DataAbstract* temp=m_data->deepCopy();
414     return new Data(temp);
415     }
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     #ifdef PASO_MPI
1127 jfenwick 1873 AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals_local,*this);
1128 ksteube 1312 // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1129     double *tmp = new double[dataPointSize];
1130     double *tmp_local = new double[dataPointSize];
1131     for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1132     MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1133     for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1134     delete[] tmp;
1135     delete[] tmp_local;
1136     #else
1137 jfenwick 1872 AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);
1138 ksteube 1312 #endif
1139 jgs 94
1140     //
1141     // create the numeric array to be returned
1142     // and load the array with the integral values
1143     boost::python::numeric::array bp_array(1.0);
1144     if (rank==0) {
1145 jgs 108 bp_array.resize(1);
1146 jgs 94 index = 0;
1147     bp_array[0] = integrals[index];
1148     }
1149     if (rank==1) {
1150     bp_array.resize(shape[0]);
1151     for (int i=0; i<shape[0]; i++) {
1152     index = i;
1153     bp_array[i] = integrals[index];
1154     }
1155     }
1156     if (rank==2) {
1157 gross 436 bp_array.resize(shape[0],shape[1]);
1158     for (int i=0; i<shape[0]; i++) {
1159     for (int j=0; j<shape[1]; j++) {
1160     index = i + shape[0] * j;
1161     bp_array[make_tuple(i,j)] = integrals[index];
1162     }
1163     }
1164 jgs 94 }
1165     if (rank==3) {
1166     bp_array.resize(shape[0],shape[1],shape[2]);
1167     for (int i=0; i<shape[0]; i++) {
1168     for (int j=0; j<shape[1]; j++) {
1169     for (int k=0; k<shape[2]; k++) {
1170     index = i + shape[0] * ( j + shape[1] * k );
1171 gross 436 bp_array[make_tuple(i,j,k)] = integrals[index];
1172 jgs 94 }
1173     }
1174     }
1175     }
1176     if (rank==4) {
1177     bp_array.resize(shape[0],shape[1],shape[2],shape[3]);
1178     for (int i=0; i<shape[0]; i++) {
1179     for (int j=0; j<shape[1]; j++) {
1180     for (int k=0; k<shape[2]; k++) {
1181     for (int l=0; l<shape[3]; l++) {
1182     index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
1183 gross 436 bp_array[make_tuple(i,j,k,l)] = integrals[index];
1184 jgs 94 }
1185     }
1186     }
1187     }
1188     }
1189    
1190     //
1191     // return the loaded array
1192     return bp_array;
1193     }
1194    
1195     Data
1196     Data::sin() const
1197     {
1198 jfenwick 2005 if (isLazy())
1199     {
1200     DataLazy* c=new DataLazy(borrowDataPtr(),SIN);
1201     return Data(c);
1202     }
1203 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1204 jgs 94 }
1205    
1206     Data
1207     Data::cos() const
1208     {
1209 jfenwick 2005 if (isLazy())
1210     {
1211     DataLazy* c=new DataLazy(borrowDataPtr(),COS);
1212     return Data(c);
1213     }
1214 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1215 jgs 94 }
1216    
1217     Data
1218     Data::tan() const
1219     {
1220 jfenwick 2005 if (isLazy())
1221     {
1222     DataLazy* c=new DataLazy(borrowDataPtr(),TAN);
1223     return Data(c);
1224     }
1225 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1226 jgs 94 }
1227    
1228     Data
1229 jgs 150 Data::asin() const
1230     {
1231 jfenwick 2005 if (isLazy())
1232     {
1233     DataLazy* c=new DataLazy(borrowDataPtr(),ASIN);
1234     return Data(c);
1235     }
1236 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1237 jgs 150 }
1238    
1239     Data
1240     Data::acos() const
1241     {
1242 jfenwick 2005 if (isLazy())
1243     {
1244     DataLazy* c=new DataLazy(borrowDataPtr(),ACOS);
1245     return Data(c);
1246     }
1247 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1248 jgs 150 }
1249    
1250 phornby 1026
1251 jgs 150 Data
1252     Data::atan() const
1253     {
1254 jfenwick 2005 if (isLazy())
1255     {
1256     DataLazy* c=new DataLazy(borrowDataPtr(),ATAN);
1257     return Data(c);
1258     }
1259 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1260 jgs 150 }
1261    
1262     Data
1263     Data::sinh() const
1264     {
1265 jfenwick 2005 if (isLazy())
1266     {
1267     DataLazy* c=new DataLazy(borrowDataPtr(),SINH);
1268     return Data(c);
1269     }
1270     return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1271 jgs 150 }
1272    
1273     Data
1274     Data::cosh() const
1275     {
1276 jfenwick 2005 if (isLazy())
1277     {
1278     DataLazy* c=new DataLazy(borrowDataPtr(),COSH);
1279     return Data(c);
1280     }
1281     return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1282 jgs 150 }
1283    
1284     Data
1285     Data::tanh() const
1286     {
1287 jfenwick 2005 if (isLazy())
1288     {
1289     DataLazy* c=new DataLazy(borrowDataPtr(),TANH);
1290     return Data(c);
1291     }
1292     return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1293 jgs 150 }
1294    
1295 phornby 1026
1296 jgs 150 Data
1297 phornby 1026 Data::erf() const
1298     {
1299 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1300 gross 1028 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1301     #else
1302 jfenwick 2005 if (isLazy())
1303     {
1304     DataLazy* c=new DataLazy(borrowDataPtr(),ERF);
1305     return Data(c);
1306     }
1307 matt 1334 return C_TensorUnaryOperation(*this, ::erf);
1308 phornby 1026 #endif
1309     }
1310    
1311     Data
1312 jgs 150 Data::asinh() const
1313     {
1314 jfenwick 2005 if (isLazy())
1315     {
1316     DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);
1317     return Data(c);
1318     }
1319 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1320 matt 1334 return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1321 phornby 1032 #else
1322 matt 1334 return C_TensorUnaryOperation(*this, ::asinh);
1323 phornby 1032 #endif
1324 jgs 150 }
1325    
1326     Data
1327     Data::acosh() const
1328     {
1329 jfenwick 2005 if (isLazy())
1330     {
1331     DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);
1332     return Data(c);
1333     }
1334 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1335 matt 1334 return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1336 phornby 1032 #else
1337 matt 1334 return C_TensorUnaryOperation(*this, ::acosh);
1338 phornby 1032 #endif
1339 jgs 150 }
1340    
1341     Data
1342     Data::atanh() const
1343     {
1344 jfenwick 2005 if (isLazy())
1345     {
1346     DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);
1347     return Data(c);
1348     }
1349 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1350 matt 1334 return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1351 phornby 1032 #else
1352 matt 1334 return C_TensorUnaryOperation(*this, ::atanh);
1353 phornby 1032 #endif
1354 jgs 150 }
1355    
1356     Data
1357 gross 286 Data::log10() const
1358 jfenwick 2005 { if (isLazy())
1359     {
1360     DataLazy* c=new DataLazy(borrowDataPtr(),LOG10);
1361     return Data(c);
1362     }
1363 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1364 jgs 94 }
1365    
1366     Data
1367 gross 286 Data::log() const
1368 jgs 94 {
1369 jfenwick 2005 if (isLazy())
1370     {
1371     DataLazy* c=new DataLazy(borrowDataPtr(),LOG);
1372     return Data(c);
1373     }
1374 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1375 jgs 94 }
1376    
1377 jgs 106 Data
1378     Data::sign() const
1379 jgs 94 {
1380 jfenwick 2005 if (isLazy())
1381     {
1382     DataLazy* c=new DataLazy(borrowDataPtr(),SIGN);
1383     return Data(c);
1384     }
1385 matt 1334 return C_TensorUnaryOperation(*this, escript::fsign);
1386 jgs 94 }
1387    
1388 jgs 106 Data
1389     Data::abs() const
1390 jgs 94 {
1391 jfenwick 2005 if (isLazy())
1392     {
1393     DataLazy* c=new DataLazy(borrowDataPtr(),ABS);
1394     return Data(c);
1395     }
1396 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1397 jgs 94 }
1398    
1399 jgs 106 Data
1400     Data::neg() const
1401 jgs 94 {
1402 jfenwick 2005 if (isLazy())
1403     {
1404     DataLazy* c=new DataLazy(borrowDataPtr(),NEG);
1405     return Data(c);
1406     }
1407 matt 1334 return C_TensorUnaryOperation(*this, negate<double>());
1408 jgs 94 }
1409    
1410 jgs 102 Data
1411 jgs 106 Data::pos() const
1412 jgs 94 {
1413 jfenwick 2005 // not doing lazy check here is deliberate.
1414     // since a deep copy of lazy data should be cheap, I'll just let it happen now
1415 jgs 148 Data result;
1416     // perform a deep copy
1417     result.copy(*this);
1418     return result;
1419 jgs 102 }
1420    
1421     Data
1422 jgs 106 Data::exp() const
1423 jfenwick 2005 {
1424     if (isLazy())
1425     {
1426     DataLazy* c=new DataLazy(borrowDataPtr(),EXP);
1427     return Data(c);
1428     }
1429 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1430 jgs 102 }
1431    
1432     Data
1433 jgs 106 Data::sqrt() const
1434 jgs 102 {
1435 jfenwick 2005 if (isLazy())
1436     {
1437     DataLazy* c=new DataLazy(borrowDataPtr(),SQRT);
1438     return Data(c);
1439     }
1440 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1441 jgs 102 }
1442    
1443 jgs 106 double
1444 jfenwick 2005 Data::Lsup_const() const
1445 jgs 102 {
1446 jfenwick 2005 if (isLazy())
1447     {
1448     throw DataException("Error - cannot compute Lsup for constant lazy data.");
1449     }
1450     return LsupWorker();
1451     }
1452    
1453     double
1454     Data::Lsup()
1455     {
1456     if (isLazy())
1457     {
1458 jfenwick 2085 resolve();
1459 jfenwick 2005 }
1460     return LsupWorker();
1461     }
1462    
1463     double
1464     Data::sup_const() const
1465     {
1466     if (isLazy())
1467     {
1468     throw DataException("Error - cannot compute sup for constant lazy data.");
1469     }
1470     return supWorker();
1471     }
1472    
1473     double
1474     Data::sup()
1475     {
1476     if (isLazy())
1477     {
1478 jfenwick 2085 resolve();
1479 jfenwick 2005 }
1480     return supWorker();
1481     }
1482    
1483     double
1484     Data::inf_const() const
1485     {
1486     if (isLazy())
1487     {
1488     throw DataException("Error - cannot compute inf for constant lazy data.");
1489     }
1490     return infWorker();
1491     }
1492    
1493     double
1494     Data::inf()
1495     {
1496     if (isLazy())
1497     {
1498 jfenwick 2085 resolve();
1499 jfenwick 2005 }
1500     return infWorker();
1501     }
1502    
1503     double
1504     Data::LsupWorker() const
1505     {
1506 phornby 1455 double localValue;
1507 jgs 106 //
1508     // set the initial absolute maximum value to zero
1509 bcumming 751
1510 jgs 147 AbsMax abs_max_func;
1511 bcumming 751 localValue = algorithm(abs_max_func,0);
1512     #ifdef PASO_MPI
1513 phornby 1455 double globalValue;
1514 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1515     return globalValue;
1516     #else
1517     return localValue;
1518     #endif
1519 jgs 117 }
1520    
1521     double
1522 jfenwick 2005 Data::supWorker() const
1523 jgs 102 {
1524 phornby 1455 double localValue;
1525 jgs 106 //
1526     // set the initial maximum value to min possible double
1527 jgs 147 FMax fmax_func;
1528 bcumming 751 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1529     #ifdef PASO_MPI
1530 phornby 1455 double globalValue;
1531 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1532     return globalValue;
1533     #else
1534     return localValue;
1535     #endif
1536 jgs 102 }
1537    
1538 jgs 106 double
1539 jfenwick 2005 Data::infWorker() const
1540 jgs 102 {
1541 phornby 1455 double localValue;
1542 jgs 106 //
1543     // set the initial minimum value to max possible double
1544 jgs 147 FMin fmin_func;
1545 bcumming 751 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1546     #ifdef PASO_MPI
1547 phornby 1455 double globalValue;
1548 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1549     return globalValue;
1550     #else
1551     return localValue;
1552     #endif
1553 jgs 102 }
1554    
1555 bcumming 751 /* TODO */
1556     /* global reduction */
1557 jgs 102 Data
1558 jgs 106 Data::maxval() const
1559 jgs 102 {
1560 jfenwick 2037 if (isLazy())
1561     {
1562     Data temp(*this); // to get around the fact that you can't resolve a const Data
1563     temp.resolve();
1564     return temp.maxval();
1565     }
1566 jgs 113 //
1567     // set the initial maximum value to min possible double
1568 jgs 147 FMax fmax_func;
1569     return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1570 jgs 102 }
1571    
1572     Data
1573 jgs 106 Data::minval() const
1574 jgs 102 {
1575 jfenwick 2037 if (isLazy())
1576     {
1577     Data temp(*this); // to get around the fact that you can't resolve a const Data
1578     temp.resolve();
1579     return temp.minval();
1580     }
1581 jgs 113 //
1582     // set the initial minimum value to max possible double
1583 jgs 147 FMin fmin_func;
1584     return dp_algorithm(fmin_func,numeric_limits<double>::max());
1585 jgs 102 }
1586    
1587 jgs 123 Data
1588 gross 804 Data::swapaxes(const int axis0, const int axis1) const
1589 jgs 123 {
1590 jfenwick 2085 if (isLazy())
1591     {
1592     Data temp(*this);
1593     temp.resolve();
1594     return temp.swapaxes(axis0,axis1);
1595     }
1596 gross 804 int axis0_tmp,axis1_tmp;
1597 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1598     DataTypes::ShapeType ev_shape;
1599 gross 800 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1600     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1601     int rank=getDataPointRank();
1602 gross 804 if (rank<2) {
1603     throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1604 gross 800 }
1605 gross 804 if (axis0<0 || axis0>rank-1) {
1606     throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1607     }
1608     if (axis1<0 || axis1>rank-1) {
1609     throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1610     }
1611     if (axis0 == axis1) {
1612     throw DataException("Error - Data::swapaxes: axis indices must be different.");
1613     }
1614     if (axis0 > axis1) {
1615     axis0_tmp=axis1;
1616     axis1_tmp=axis0;
1617     } else {
1618     axis0_tmp=axis0;
1619     axis1_tmp=axis1;
1620     }
1621 gross 800 for (int i=0; i<rank; i++) {
1622 gross 804 if (i == axis0_tmp) {
1623 ksteube 1312 ev_shape.push_back(s[axis1_tmp]);
1624 gross 804 } else if (i == axis1_tmp) {
1625 ksteube 1312 ev_shape.push_back(s[axis0_tmp]);
1626 gross 800 } else {
1627 ksteube 1312 ev_shape.push_back(s[i]);
1628 gross 800 }
1629     }
1630     Data ev(0.,ev_shape,getFunctionSpace());
1631     ev.typeMatchRight(*this);
1632 gross 804 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1633 gross 800 return ev;
1634    
1635 jgs 123 }
1636    
1637     Data
1638 ksteube 775 Data::symmetric() const
1639 jgs 123 {
1640 ksteube 775 // check input
1641 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1642 ksteube 775 if (getDataPointRank()==2) {
1643 ksteube 1312 if(s[0] != s[1])
1644 ksteube 775 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1645     }
1646     else if (getDataPointRank()==4) {
1647     if(!(s[0] == s[2] && s[1] == s[3]))
1648     throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1649     }
1650     else {
1651     throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1652     }
1653 jfenwick 2005 if (isLazy())
1654     {
1655 jfenwick 2037 DataLazy* c=new DataLazy(borrowDataPtr(),SYM);
1656     return Data(c);
1657 jfenwick 2005 }
1658 ksteube 775 Data ev(0.,getDataPointShape(),getFunctionSpace());
1659     ev.typeMatchRight(*this);
1660     m_data->symmetric(ev.m_data.get());
1661     return ev;
1662     }
1663    
1664     Data
1665     Data::nonsymmetric() const
1666     {
1667 jfenwick 2005 if (isLazy())
1668     {
1669 jfenwick 2037 DataLazy* c=new DataLazy(borrowDataPtr(),NSYM);
1670     return Data(c);
1671 jfenwick 2005 }
1672 ksteube 775 // check input
1673 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1674 ksteube 775 if (getDataPointRank()==2) {
1675 ksteube 1312 if(s[0] != s[1])
1676 ksteube 775 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1677 jfenwick 1796 DataTypes::ShapeType ev_shape;
1678 ksteube 775 ev_shape.push_back(s[0]);
1679     ev_shape.push_back(s[1]);
1680     Data ev(0.,ev_shape,getFunctionSpace());
1681     ev.typeMatchRight(*this);
1682     m_data->nonsymmetric(ev.m_data.get());
1683     return ev;
1684     }
1685     else if (getDataPointRank()==4) {
1686     if(!(s[0] == s[2] && s[1] == s[3]))
1687     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1688 jfenwick 1796 DataTypes::ShapeType ev_shape;
1689 ksteube 775 ev_shape.push_back(s[0]);
1690     ev_shape.push_back(s[1]);
1691     ev_shape.push_back(s[2]);
1692     ev_shape.push_back(s[3]);
1693     Data ev(0.,ev_shape,getFunctionSpace());
1694     ev.typeMatchRight(*this);
1695     m_data->nonsymmetric(ev.m_data.get());
1696     return ev;
1697     }
1698     else {
1699     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1700     }
1701     }
1702    
1703     Data
1704 gross 800 Data::trace(int axis_offset) const
1705 jfenwick 2084 {
1706 jfenwick 2005 if (isLazy())
1707     {
1708 jfenwick 2084 DataLazy* c=new DataLazy(borrowDataPtr(),TRACE,axis_offset);
1709     return Data(c);
1710 jfenwick 2005 }
1711 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1712 ksteube 775 if (getDataPointRank()==2) {
1713 jfenwick 1796 DataTypes::ShapeType ev_shape;
1714 ksteube 775 Data ev(0.,ev_shape,getFunctionSpace());
1715     ev.typeMatchRight(*this);
1716 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1717 ksteube 775 return ev;
1718     }
1719     if (getDataPointRank()==3) {
1720 jfenwick 1796 DataTypes::ShapeType ev_shape;
1721 ksteube 775 if (axis_offset==0) {
1722     int s2=s[2];
1723     ev_shape.push_back(s2);
1724     }
1725     else if (axis_offset==1) {
1726     int s0=s[0];
1727     ev_shape.push_back(s0);
1728     }
1729     Data ev(0.,ev_shape,getFunctionSpace());
1730     ev.typeMatchRight(*this);
1731 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1732 ksteube 775 return ev;
1733     }
1734     if (getDataPointRank()==4) {
1735 jfenwick 1796 DataTypes::ShapeType ev_shape;
1736 ksteube 775 if (axis_offset==0) {
1737     ev_shape.push_back(s[2]);
1738     ev_shape.push_back(s[3]);
1739     }
1740     else if (axis_offset==1) {
1741     ev_shape.push_back(s[0]);
1742     ev_shape.push_back(s[3]);
1743     }
1744     else if (axis_offset==2) {
1745     ev_shape.push_back(s[0]);
1746     ev_shape.push_back(s[1]);
1747     }
1748     Data ev(0.,ev_shape,getFunctionSpace());
1749     ev.typeMatchRight(*this);
1750 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1751 ksteube 775 return ev;
1752     }
1753     else {
1754 gross 800 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1755 ksteube 775 }
1756     }
1757    
1758     Data
1759     Data::transpose(int axis_offset) const
1760 jfenwick 2005 {
1761     if (isLazy())
1762     {
1763 jfenwick 2084 DataLazy* c=new DataLazy(borrowDataPtr(),TRANS,axis_offset);
1764     return Data(c);
1765 jfenwick 2005 }
1766 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1767     DataTypes::ShapeType ev_shape;
1768 ksteube 775 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1769     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1770     int rank=getDataPointRank();
1771     if (axis_offset<0 || axis_offset>rank) {
1772     throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1773     }
1774     for (int i=0; i<rank; i++) {
1775     int index = (axis_offset+i)%rank;
1776     ev_shape.push_back(s[index]); // Append to new shape
1777     }
1778     Data ev(0.,ev_shape,getFunctionSpace());
1779     ev.typeMatchRight(*this);
1780     m_data->transpose(ev.m_data.get(), axis_offset);
1781     return ev;
1782 jgs 123 }
1783    
1784 gross 576 Data
1785     Data::eigenvalues() const
1786     {
1787 jfenwick 2005 if (isLazy())
1788     {
1789     Data temp(*this); // to get around the fact that you can't resolve a const Data
1790     temp.resolve();
1791     return temp.eigenvalues();
1792     }
1793 gross 576 // check input
1794 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1795 ksteube 1312 if (getDataPointRank()!=2)
1796 gross 576 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1797 ksteube 1312 if(s[0] != s[1])
1798 gross 576 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1799     // create return
1800 jfenwick 1796 DataTypes::ShapeType ev_shape(1,s[0]);
1801 gross 576 Data ev(0.,ev_shape,getFunctionSpace());
1802     ev.typeMatchRight(*this);
1803     m_data->eigenvalues(ev.m_data.get());
1804     return ev;
1805     }
1806    
1807 jgs 121 const boost::python::tuple
1808 gross 576 Data::eigenvalues_and_eigenvectors(const double tol) const
1809     {
1810 jfenwick 2005 if (isLazy())
1811     {
1812     Data temp(*this); // to get around the fact that you can't resolve a const Data
1813     temp.resolve();
1814     return temp.eigenvalues_and_eigenvectors(tol);
1815     }
1816 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1817 ksteube 1312 if (getDataPointRank()!=2)
1818 gross 576 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1819 ksteube 1312 if(s[0] != s[1])
1820 gross 576 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1821     // create return
1822 jfenwick 1796 DataTypes::ShapeType ev_shape(1,s[0]);
1823 gross 576 Data ev(0.,ev_shape,getFunctionSpace());
1824     ev.typeMatchRight(*this);
1825 jfenwick 1796 DataTypes::ShapeType V_shape(2,s[0]);
1826 gross 576 Data V(0.,V_shape,getFunctionSpace());
1827     V.typeMatchRight(*this);
1828     m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1829     return make_tuple(boost::python::object(ev),boost::python::object(V));
1830     }
1831    
1832     const boost::python::tuple
1833 gross 921 Data::minGlobalDataPoint() const
1834 jgs 121 {
1835 gross 921 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1836 jgs 148 // abort (for unknown reasons) if there are openmp directives with it in the
1837     // surrounding function
1838    
1839     int DataPointNo;
1840 gross 921 int ProcNo;
1841     calc_minGlobalDataPoint(ProcNo,DataPointNo);
1842     return make_tuple(ProcNo,DataPointNo);
1843 jgs 148 }
1844    
1845     void
1846 gross 921 Data::calc_minGlobalDataPoint(int& ProcNo,
1847     int& DataPointNo) const
1848 jgs 148 {
1849 jfenwick 2005 if (isLazy())
1850     {
1851     Data temp(*this); // to get around the fact that you can't resolve a const Data
1852     temp.resolve();
1853     return temp.calc_minGlobalDataPoint(ProcNo,DataPointNo);
1854     }
1855 jgs 148 int i,j;
1856     int lowi=0,lowj=0;
1857     double min=numeric_limits<double>::max();
1858    
1859 jgs 121 Data temp=minval();
1860    
1861     int numSamples=temp.getNumSamples();
1862     int numDPPSample=temp.getNumDataPointsPerSample();
1863    
1864 jgs 148 double next,local_min;
1865 jfenwick 1946 int local_lowi=0,local_lowj=0;
1866 jgs 121
1867 caltinay 2081 #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
1868 jgs 148 {
1869     local_min=min;
1870     #pragma omp for private(i,j) schedule(static)
1871     for (i=0; i<numSamples; i++) {
1872     for (j=0; j<numDPPSample; j++) {
1873 jfenwick 1796 next=temp.getDataAtOffset(temp.getDataOffset(i,j));
1874 jgs 148 if (next<local_min) {
1875     local_min=next;
1876     local_lowi=i;
1877     local_lowj=j;
1878     }
1879 jgs 121 }
1880     }
1881 jgs 148 #pragma omp critical
1882     if (local_min<min) {
1883     min=local_min;
1884     lowi=local_lowi;
1885     lowj=local_lowj;
1886     }
1887 jgs 121 }
1888    
1889 bcumming 782 #ifdef PASO_MPI
1890     // determine the processor on which the minimum occurs
1891 jfenwick 1796 next = temp.getDataPoint(lowi,lowj);
1892 bcumming 782 int lowProc = 0;
1893     double *globalMins = new double[get_MPISize()+1];
1894 caltinay 2081 int error;
1895     error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1896 ksteube 1312
1897 bcumming 782 if( get_MPIRank()==0 ){
1898     next = globalMins[lowProc];
1899     for( i=1; i<get_MPISize(); i++ )
1900     if( next>globalMins[i] ){
1901     lowProc = i;
1902     next = globalMins[i];
1903     }
1904     }
1905     MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1906    
1907     delete [] globalMins;
1908     ProcNo = lowProc;
1909 bcumming 790 #else
1910     ProcNo = 0;
1911 bcumming 782 #endif
1912 gross 921 DataPointNo = lowj + lowi * numDPPSample;
1913 jgs 121 }
1914    
1915 jgs 104 void
1916     Data::saveDX(std::string fileName) const
1917     {
1918 jfenwick 1803 if (isEmpty())
1919     {
1920     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1921     }
1922 jfenwick 2005 if (isLazy())
1923     {
1924     Data temp(*this); // to get around the fact that you can't resolve a const Data
1925     temp.resolve();
1926     temp.saveDX(fileName);
1927     return;
1928     }
1929 jgs 153 boost::python::dict args;
1930     args["data"]=boost::python::object(this);
1931 jfenwick 1872 getDomain()->saveDX(fileName,args);
1932 jgs 104 return;
1933     }
1934    
1935 jgs 110 void
1936     Data::saveVTK(std::string fileName) const
1937     {
1938 jfenwick 1803 if (isEmpty())
1939     {
1940     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1941     }
1942 jfenwick 2005 if (isLazy())
1943     {
1944     Data temp(*this); // to get around the fact that you can't resolve a const Data
1945     temp.resolve();
1946     temp.saveVTK(fileName);
1947     return;
1948     }
1949 jgs 153 boost::python::dict args;
1950     args["data"]=boost::python::object(this);
1951 jfenwick 1872 getDomain()->saveVTK(fileName,args);
1952 jgs 110 return;
1953     }
1954    
1955 jgs 102 Data&
1956     Data::operator+=(const Data& right)
1957     {
1958 gross 783 if (isProtected()) {
1959     throw DataException("Error - attempt to update protected Data object.");
1960     }
1961 jfenwick 2005 if (isLazy() || right.isLazy())
1962     {
1963     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),ADD); // for lazy + is equivalent to +=
1964     m_data=c->getPtr();
1965     return (*this);
1966     }
1967     else
1968     {
1969     binaryOp(right,plus<double>());
1970     return (*this);
1971     }
1972 jgs 94 }
1973    
1974 jgs 102 Data&
1975     Data::operator+=(const boost::python::object& right)
1976 jgs 94 {
1977 jfenwick 2005 if (isProtected()) {
1978     throw DataException("Error - attempt to update protected Data object.");
1979     }
1980 gross 854 Data tmp(right,getFunctionSpace(),false);
1981 jfenwick 2005 if (isLazy())
1982     {
1983     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),ADD); // for lazy + is equivalent to +=
1984     m_data=c->getPtr();
1985     return (*this);
1986     }
1987     else
1988     {
1989     binaryOp(tmp,plus<double>());
1990     return (*this);
1991     }
1992 jgs 94 }
1993 jfenwick 2005
1994     // Hmmm, operator= makes a deep copy but the copy constructor does not?
1995 ksteube 1312 Data&
1996     Data::operator=(const Data& other)
1997     {
1998     copy(other);
1999     return (*this);
2000     }
2001 jgs 94
2002 jgs 102 Data&
2003     Data::operator-=(const Data& right)
2004 jgs 94 {
2005 gross 783 if (isProtected()) {
2006     throw DataException("Error - attempt to update protected Data object.");
2007     }
2008 jfenwick 2005 if (isLazy() || right.isLazy())
2009     {
2010     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),SUB); // for lazy - is equivalent to -=
2011     m_data=c->getPtr();
2012     return (*this);
2013     }
2014     else
2015     {
2016     binaryOp(right,minus<double>());
2017     return (*this);
2018     }
2019 jgs 94 }
2020    
2021 jgs 102 Data&
2022     Data::operator-=(const boost::python::object& right)
2023 jgs 94 {
2024 jfenwick 2005 if (isProtected()) {
2025     throw DataException("Error - attempt to update protected Data object.");
2026     }
2027 gross 854 Data tmp(right,getFunctionSpace(),false);
2028 jfenwick 2005 if (isLazy())
2029     {
2030     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),SUB); // for lazy - is equivalent to -=
2031     m_data=c->getPtr();
2032     return (*this);
2033     }
2034     else
2035     {
2036     binaryOp(tmp,minus<double>());
2037     return (*this);
2038     }
2039 jgs 94 }
2040    
2041 jgs 102 Data&
2042     Data::operator*=(const Data& right)
2043 jgs 94 {
2044 gross 783 if (isProtected()) {
2045     throw DataException("Error - attempt to update protected Data object.");
2046     }
2047 jfenwick 2005 if (isLazy() || right.isLazy())
2048     {
2049     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),MUL); // for lazy * is equivalent to *=
2050     m_data=c->getPtr();
2051     return (*this);
2052     }
2053     else
2054     {
2055     binaryOp(right,multiplies<double>());
2056     return (*this);
2057     }
2058 jgs 94 }
2059    
2060 jgs 102 Data&
2061     Data::operator*=(const boost::python::object& right)
2062 jfenwick 2005 {
2063     if (isProtected()) {
2064     throw DataException("Error - attempt to update protected Data object.");
2065     }
2066 gross 854 Data tmp(right,getFunctionSpace(),false);
2067 jfenwick 2005 if (isLazy())
2068     {
2069     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),MUL); // for lazy * is equivalent to *=
2070     m_data=c->getPtr();
2071     return (*this);
2072     }
2073     else
2074     {
2075     binaryOp(tmp,multiplies<double>());
2076     return (*this);
2077     }
2078 jgs 94 }
2079    
2080 jgs 102 Data&
2081     Data::operator/=(const Data& right)
2082 jgs 94 {
2083 gross 783 if (isProtected()) {
2084     throw DataException("Error - attempt to update protected Data object.");
2085     }
2086 jfenwick 2005 if (isLazy() || right.isLazy())
2087     {
2088     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),DIV); // for lazy / is equivalent to /=
2089     m_data=c->getPtr();
2090     return (*this);
2091     }
2092     else
2093     {
2094     binaryOp(right,divides<double>());
2095     return (*this);
2096     }
2097 jgs 94 }
2098    
2099 jgs 102 Data&
2100     Data::operator/=(const boost::python::object& right)
2101 jgs 94 {
2102 jfenwick 2005 if (isProtected()) {
2103     throw DataException("Error - attempt to update protected Data object.");
2104     }
2105 gross 854 Data tmp(right,getFunctionSpace(),false);
2106 jfenwick 2005 if (isLazy())
2107     {
2108     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),DIV); // for lazy / is equivalent to /=
2109     m_data=c->getPtr();
2110     return (*this);
2111     }
2112     else
2113     {
2114     binaryOp(tmp,divides<double>());
2115     return (*this);
2116     }
2117 jgs 94 }
2118    
2119 jgs 102 Data
2120 gross 699 Data::rpowO(const boost::python::object& left) const
2121     {
2122     Data left_d(left,*this);
2123     return left_d.powD(*this);
2124     }
2125    
2126     Data
2127 jgs 102 Data::powO(const boost::python::object& right) const
2128 jgs 94 {
2129 gross 854 Data tmp(right,getFunctionSpace(),false);
2130     return powD(tmp);
2131 jgs 94 }
2132    
2133 jgs 102 Data
2134     Data::powD(const Data& right) const
2135 jgs 94 {
2136 jfenwick 2005 if (isLazy() || right.isLazy())
2137     {
2138     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),POW);
2139     return Data(c);
2140     }
2141 matt 1350 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2142 jgs 94 }
2143    
2144     //
2145 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2146 jgs 102 Data
2147     escript::operator+(const Data& left, const Data& right)
2148 jgs 94 {
2149 jfenwick 2005 if (left.isLazy() || right.isLazy())
2150     {
2151     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);
2152     return Data(c);
2153     }
2154 matt 1327 return C_TensorBinaryOperation(left, right, plus<double>());
2155 jgs 94 }
2156    
2157     //
2158 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2159 jgs 102 Data
2160     escript::operator-(const Data& left, const Data& right)
2161 jgs 94 {
2162 jfenwick 2005 if (left.isLazy() || right.isLazy())
2163     {
2164     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);
2165     return Data(c);
2166     }
2167 matt 1327 return C_TensorBinaryOperation(left, right, minus<double>());
2168 jgs 94 }
2169    
2170     //
2171 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2172 jgs 102 Data
2173     escript::operator*(const Data& left, const Data& right)
2174 jgs 94 {
2175 jfenwick 2005 if (left.isLazy() || right.isLazy())
2176     {
2177     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);
2178     return Data(c);
2179     }
2180 matt 1327 return C_TensorBinaryOperation(left, right, multiplies<double>());
2181 jgs 94 }
2182    
2183     //
2184 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2185 jgs 102 Data
2186     escript::operator/(const Data& left, const Data& right)
2187 jgs 94 {
2188 jfenwick 2005 if (left.isLazy() || right.isLazy())
2189     {
2190     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),DIV);
2191     return Data(c);
2192     }
2193 matt 1327 return C_TensorBinaryOperation(left, right, divides<double>());
2194 jgs 94 }
2195    
2196     //
2197 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2198 jgs 102 Data
2199     escript::operator+(const Data& left, const boost::python::object& right)
2200 jgs 94 {
2201 jfenwick 2005 if (left.isLazy())
2202     {
2203     DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),ADD);
2204     return Data(c);
2205     }
2206 gross 854 return left+Data(right,left.getFunctionSpace(),false);
2207 jgs 94 }
2208    
2209     //
2210 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2211 jgs 102 Data
2212     escript::operator-(const Data& left, const boost::python::object& right)
2213 jgs 94 {
2214 jfenwick 2005 if (left.isLazy())
2215     {
2216     DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),SUB);
2217     return Data(c);
2218     }
2219 gross 854 return left-Data(right,left.getFunctionSpace(),false);
2220 jgs 94 }
2221    
2222     //
2223 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2224 jgs 102 Data
2225     escript::operator*(const Data& left, const boost::python::object& right)
2226 jgs 94 {
2227 jfenwick 2005 if (left.isLazy())
2228     {
2229     DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),MUL);
2230     return Data(c);
2231     }
2232 gross 854 return left*Data(right,left.getFunctionSpace(),false);
2233 jgs 94 }
2234    
2235     //
2236 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2237 jgs 102 Data
2238     escript::operator/(const Data& left, const boost::python::object& right)
2239 jgs 94 {
2240 jfenwick 2005 if (left.isLazy())
2241     {
2242     DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),DIV);
2243     return Data(c);
2244     }
2245 gross 854 return left/Data(right,left.getFunctionSpace(),false);
2246 jgs 94 }
2247    
2248     //
2249 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2250 jgs 102 Data
2251     escript::operator+(const boost::python::object& left, const Data& right)
2252 jgs 94 {
2253 jfenwick 2005 if (right.isLazy())
2254     {
2255     DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),ADD);
2256     return Data(c);
2257     }
2258 gross 854 return Data(left,right.getFunctionSpace(),false)+right;
2259 jgs 94 }
2260    
2261     //
2262 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2263 jgs 102 Data
2264     escript::operator-(const boost::python::object& left, const Data& right)
2265 jgs 94 {
2266 jfenwick 2005 if (right.isLazy())
2267     {
2268     DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),SUB);
2269     return Data(c);
2270     }
2271 gross 854 return Data(left,right.getFunctionSpace(),false)-right;
2272 jgs 94 }
2273    
2274     //
2275 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2276 jgs 102 Data
2277     escript::operator*(const boost::python::object& left, const Data& right)
2278 jgs 94 {
2279 jfenwick 2005 if (right.isLazy())
2280     {
2281     DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),MUL);
2282     return Data(c);
2283     }
2284 gross 854 return Data(left,right.getFunctionSpace(),false)*right;
2285 jgs 94 }
2286    
2287     //
2288 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2289 jgs 102 Data
2290     escript::operator/(const boost::python::object& left, const Data& right)
2291 jgs 94 {
2292 jfenwick 2005 if (right.isLazy())
2293     {
2294     DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),DIV);
2295     return Data(c);
2296     }
2297 gross 854 return Data(left,right.getFunctionSpace(),false)/right;
2298 jgs 94 }
2299    
2300 jgs 102
2301 bcumming 751 /* TODO */
2302     /* global reduction */
2303 jgs 102 Data
2304 ksteube 1312 Data::getItem(const boost::python::object& key) const
2305 jgs 94 {
2306    
2307 jfenwick 1796 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2308 jgs 94
2309 jfenwick 1796 if (slice_region.size()!=getDataPointRank()) {
2310 jgs 102 throw DataException("Error - slice size does not match Data rank.");
2311 jgs 94 }
2312    
2313 jgs 102 return getSlice(slice_region);
2314 jgs 94 }
2315    
2316 bcumming 751 /* TODO */
2317     /* global reduction */
2318 jgs 94 Data
2319 jfenwick 1796 Data::getSlice(const DataTypes::RegionType& region) const
2320 jgs 94 {
2321 jgs 102 return Data(*this,region);
2322 jgs 94 }
2323    
2324 bcumming 751 /* TODO */
2325     /* global reduction */
2326 jgs 94 void
2327 jgs 102 Data::setItemO(const boost::python::object& key,
2328     const boost::python::object& value)
2329 jgs 94 {
2330 jgs 102 Data tempData(value,getFunctionSpace());
2331     setItemD(key,tempData);
2332     }
2333    
2334     void
2335     Data::setItemD(const boost::python::object& key,
2336     const Data& value)
2337     {
2338 jfenwick 1796 // const DataArrayView& view=getPointDataView();
2339 jgs 123
2340 jfenwick 1796 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2341     if (slice_region.size()!=getDataPointRank()) {
2342 jgs 94 throw DataException("Error - slice size does not match Data rank.");
2343     }
2344 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
2345     setSlice(Data(value,getFunctionSpace()),slice_region);
2346     } else {
2347     setSlice(value,slice_region);
2348     }
2349 jgs 94 }
2350    
2351     void
2352 jgs 102 Data::setSlice(const Data& value,
2353 jfenwick 1796 const DataTypes::RegionType& region)
2354 jgs 94 {
2355 gross 783 if (isProtected()) {
2356     throw DataException("Error - attempt to update protected Data object.");
2357     }
2358 jfenwick 2005 FORCERESOLVE;
2359     /* if (isLazy())
2360     {
2361     throw DataException("Error - setSlice not permitted on lazy data.");
2362     }*/
2363 jgs 102 Data tempValue(value);
2364     typeMatchLeft(tempValue);
2365     typeMatchRight(tempValue);
2366 jfenwick 2005 getReady()->setSlice(tempValue.m_data.get(),region);
2367 jgs 102 }
2368    
2369     void
2370     Data::typeMatchLeft(Data& right) const
2371     {
2372 jfenwick 2005 if (right.isLazy() && !isLazy())
2373     {
2374     right.resolve();
2375     }
2376 jgs 102 if (isExpanded()){
2377     right.expand();
2378     } else if (isTagged()) {
2379     if (right.isConstant()) {
2380     right.tag();
2381     }
2382     }
2383     }
2384    
2385     void
2386     Data::typeMatchRight(const Data& right)
2387     {
2388 jfenwick 2005 if (isLazy() && !right.isLazy())
2389     {
2390     resolve();
2391     }
2392 jgs 94 if (isTagged()) {
2393     if (right.isExpanded()) {
2394     expand();
2395     }
2396     } else if (isConstant()) {
2397     if (right.isExpanded()) {
2398     expand();
2399     } else if (right.isTagged()) {
2400     tag();
2401     }
2402     }
2403     }
2404    
2405     void
2406 gross 1044 Data::setTaggedValueByName(std::string name,
2407 ksteube 1312 const boost::python::object& value)
2408 gross 1044 {
2409 jfenwick 1872 if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2410 jfenwick 2005 FORCERESOLVE;
2411 jfenwick 1872 int tagKey=getFunctionSpace().getDomain()->getTag(name);
2412 gross 1044 setTaggedValue(tagKey,value);
2413     }
2414     }
2415     void
2416 jgs 94 Data::setTaggedValue(int tagKey,
2417     const boost::python::object& value)
2418     {
2419 gross 783 if (isProtected()) {
2420     throw DataException("Error - attempt to update protected Data object.");
2421     }
2422 jgs 94 //
2423     // Ensure underlying data object is of type DataTagged
2424 jfenwick 2005 FORCERESOLVE;
2425 gross 1358 if (isConstant()) tag();
2426 matt 1319 numeric::array asNumArray(value);
2427 jgs 94
2428 matt 1319 // extract the shape of the numarray
2429 jfenwick 1796 DataTypes::ShapeType tempShape;
2430 matt 1319 for (int i=0; i < asNumArray.getrank(); i++) {
2431     tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
2432     }
2433    
2434 jfenwick 1796 DataVector temp_data2;
2435     temp_data2.copyFromNumArray(asNumArray);
2436    
2437 jfenwick 2005 m_data->setTaggedValue(tagKey,tempShape, temp_data2);
2438 jgs 94 }
2439    
2440 jfenwick 1796
2441 jgs 110 void
2442 jgs 121 Data::setTaggedValueFromCPP(int tagKey,
2443 jfenwick 1796 const DataTypes::ShapeType& pointshape,
2444     const DataTypes::ValueType& value,
2445     int dataOffset)
2446 jgs 121 {
2447 gross 783 if (isProtected()) {
2448     throw DataException("Error - attempt to update protected Data object.");
2449     }
2450 jgs 121 //
2451     // Ensure underlying data object is of type DataTagged
2452 jfenwick 2005 FORCERESOLVE;
2453 gross 1358 if (isConstant()) tag();
2454 jgs 121 //
2455