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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1026 - (hide annotations)
Tue Mar 13 08:25:00 2007 UTC (12 years, 6 months ago) by phornby
File size: 82020 byte(s)
Data.cpp - indefed out the erf & inv. hyperbolics on windows.
DataAlgorithmAdapterTestCase & DataTestCase - Fix the ifndef on _WIN32 & _INTEL_COMPILER
SConstruct - chamges to the PATH so windows can find DLLs.


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