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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 854 - (hide annotations)
Thu Sep 21 05:29:42 2006 UTC (13 years, 6 months ago) by gross
File size: 82904 byte(s)
Some modifications to the binary operations +,-,*/, pow. 
The code is a bit simpler now and more efficient has there is
no reseising required now. the resizing method has been removed as
it is very, very inefficient. Even serial code should be faster now.
It is now forbidden to do an inplace update of scalar data object with an object 
of rank >0 as this is very slow (and does not make much sense). 


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