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

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

Parent Directory Parent Directory | Revision Log Revision Log


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