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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2048 - (hide annotations)
Mon Nov 17 08:46:00 2008 UTC (11 years, 4 months ago) by phornby
File size: 86143 byte(s)
Experimental commit to move the code to a windows box
on the other side of a firewall. Purpose: add support for the
imploved intelc math library on windows.


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