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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 876 - (hide annotations)
Thu Oct 19 03:50:23 2006 UTC (13 years, 1 month ago) by ksteube
File size: 83039 byte(s)
Added erf (error function) implementation

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 ksteube 876 Data::erf() const
1036     {
1037     #if defined DOPROF
1038     profData->unary++;
1039     #endif
1040     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::erf);
1041     }
1042    
1043     Data
1044 jgs 94 Data::sin() 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)::sin);
1050     }
1051    
1052     Data
1053     Data::cos() 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)::cos);
1059     }
1060    
1061     Data
1062     Data::tan() const
1063     {
1064 jgs 151 #if defined DOPROF
1065 jgs 123 profData->unary++;
1066 jgs 151 #endif
1067 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
1068     }
1069    
1070     Data
1071 jgs 150 Data::asin() 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)::asin);
1077     }
1078    
1079     Data
1080     Data::acos() 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)::acos);
1086     }
1087    
1088     Data
1089     Data::atan() 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)::atan);
1095     }
1096    
1097     Data
1098     Data::sinh() 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)::sinh);
1104     }
1105    
1106     Data
1107     Data::cosh() 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)::cosh);
1113     }
1114    
1115     Data
1116     Data::tanh() 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)::tanh);
1122     }
1123    
1124     Data
1125     Data::asinh() 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)::asinh);
1131     }
1132    
1133     Data
1134     Data::acosh() 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)::acosh);
1140     }
1141    
1142     Data
1143     Data::atanh() const
1144     {
1145 jgs 151 #if defined DOPROF
1146 jgs 150 profData->unary++;
1147 jgs 151 #endif
1148 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
1149     }
1150    
1151     Data
1152 gross 286 Data::log10() 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)::log10);
1158     }
1159    
1160     Data
1161 gross 286 Data::log() const
1162 jgs 94 {
1163 jgs 151 #if defined DOPROF
1164 jgs 123 profData->unary++;
1165 jgs 151 #endif
1166 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
1167     }
1168    
1169 jgs 106 Data
1170     Data::sign() 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,escript::fsign);
1176 jgs 94 }
1177    
1178 jgs 106 Data
1179     Data::abs() 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,(Data::UnaryDFunPtr)::fabs);
1185 jgs 94 }
1186    
1187 jgs 106 Data
1188     Data::neg() const
1189 jgs 94 {
1190 jgs 151 #if defined DOPROF
1191 jgs 123 profData->unary++;
1192 jgs 151 #endif
1193 jgs 106 return escript::unaryOp(*this,negate<double>());
1194 jgs 94 }
1195    
1196 jgs 102 Data
1197 jgs 106 Data::pos() const
1198 jgs 94 {
1199 jgs 151 #if defined DOPROF
1200 jgs 123 profData->unary++;
1201 jgs 151 #endif
1202 jgs 148 Data result;
1203     // perform a deep copy
1204     result.copy(*this);
1205     return result;
1206 jgs 102 }
1207    
1208     Data
1209 jgs 106 Data::exp() 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)::exp);
1215 jgs 102 }
1216    
1217     Data
1218 jgs 106 Data::sqrt() const
1219 jgs 102 {
1220 jgs 151 #if defined DOPROF
1221 jgs 123 profData->unary++;
1222 jgs 151 #endif
1223 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1224 jgs 102 }
1225    
1226 jgs 106 double
1227     Data::Lsup() const
1228 jgs 102 {
1229 bcumming 751 double localValue, globalValue;
1230 jgs 151 #if defined DOPROF
1231 jgs 123 profData->reduction1++;
1232 jgs 151 #endif
1233 jgs 106 //
1234     // set the initial absolute maximum value to zero
1235 bcumming 751
1236 jgs 147 AbsMax abs_max_func;
1237 bcumming 751 localValue = algorithm(abs_max_func,0);
1238     #ifdef PASO_MPI
1239     MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1240     return globalValue;
1241     #else
1242     return localValue;
1243     #endif
1244 jgs 102 }
1245    
1246 jgs 106 double
1247 jgs 117 Data::Linf() const
1248     {
1249 bcumming 751 double localValue, globalValue;
1250 jgs 151 #if defined DOPROF
1251 jgs 123 profData->reduction1++;
1252 jgs 151 #endif
1253 jgs 117 //
1254     // set the initial absolute minimum value to max double
1255 jgs 147 AbsMin abs_min_func;
1256 bcumming 751 localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1257    
1258     #ifdef PASO_MPI
1259     MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1260     return globalValue;
1261     #else
1262     return localValue;
1263     #endif
1264 jgs 117 }
1265    
1266     double
1267 jgs 106 Data::sup() const
1268 jgs 102 {
1269 bcumming 751 double localValue, globalValue;
1270 jgs 151 #if defined DOPROF
1271 jgs 123 profData->reduction1++;
1272 jgs 151 #endif
1273 jgs 106 //
1274     // set the initial maximum value to min possible double
1275 jgs 147 FMax fmax_func;
1276 bcumming 751 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1277     #ifdef PASO_MPI
1278     MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1279     return globalValue;
1280     #else
1281     return localValue;
1282     #endif
1283 jgs 102 }
1284    
1285 jgs 106 double
1286     Data::inf() const
1287 jgs 102 {
1288 bcumming 751 double localValue, globalValue;
1289 jgs 151 #if defined DOPROF
1290 jgs 123 profData->reduction1++;
1291 jgs 151 #endif
1292 jgs 106 //
1293     // set the initial minimum value to max possible double
1294 jgs 147 FMin fmin_func;
1295 bcumming 751 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1296     #ifdef PASO_MPI
1297     MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1298     return globalValue;
1299     #else
1300     return localValue;
1301     #endif
1302 jgs 102 }
1303    
1304 bcumming 751 /* TODO */
1305     /* global reduction */
1306 jgs 102 Data
1307 jgs 106 Data::maxval() const
1308 jgs 102 {
1309 jgs 151 #if defined DOPROF
1310 jgs 123 profData->reduction2++;
1311 jgs 151 #endif
1312 jgs 113 //
1313     // set the initial maximum value to min possible double
1314 jgs 147 FMax fmax_func;
1315     return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1316 jgs 102 }
1317    
1318     Data
1319 jgs 106 Data::minval() const
1320 jgs 102 {
1321 jgs 151 #if defined DOPROF
1322 jgs 123 profData->reduction2++;
1323 jgs 151 #endif
1324 jgs 113 //
1325     // set the initial minimum value to max possible double
1326 jgs 147 FMin fmin_func;
1327     return dp_algorithm(fmin_func,numeric_limits<double>::max());
1328 jgs 102 }
1329    
1330 jgs 123 Data
1331 gross 804 Data::swapaxes(const int axis0, const int axis1) const
1332 jgs 123 {
1333 gross 804 int axis0_tmp,axis1_tmp;
1334 gross 800 #if defined DOPROF
1335     profData->unary++;
1336     #endif
1337     DataArrayView::ShapeType s=getDataPointShape();
1338     DataArrayView::ShapeType ev_shape;
1339     // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1340     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1341     int rank=getDataPointRank();
1342 gross 804 if (rank<2) {
1343     throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1344 gross 800 }
1345 gross 804 if (axis0<0 || axis0>rank-1) {
1346     throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1347     }
1348     if (axis1<0 || axis1>rank-1) {
1349     throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1350     }
1351     if (axis0 == axis1) {
1352     throw DataException("Error - Data::swapaxes: axis indices must be different.");
1353     }
1354     if (axis0 > axis1) {
1355     axis0_tmp=axis1;
1356     axis1_tmp=axis0;
1357     } else {
1358     axis0_tmp=axis0;
1359     axis1_tmp=axis1;
1360     }
1361 gross 800 for (int i=0; i<rank; i++) {
1362 gross 804 if (i == axis0_tmp) {
1363     ev_shape.push_back(s[axis1_tmp]);
1364     } else if (i == axis1_tmp) {
1365     ev_shape.push_back(s[axis0_tmp]);
1366 gross 800 } else {
1367     ev_shape.push_back(s[i]);
1368     }
1369     }
1370     Data ev(0.,ev_shape,getFunctionSpace());
1371     ev.typeMatchRight(*this);
1372 gross 804 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1373 gross 800 return ev;
1374    
1375 jgs 123 }
1376    
1377     Data
1378 ksteube 775 Data::symmetric() const
1379 jgs 123 {
1380 ksteube 775 #if defined DOPROF
1381     profData->unary++;
1382     #endif
1383     // check input
1384     DataArrayView::ShapeType s=getDataPointShape();
1385     if (getDataPointRank()==2) {
1386     if(s[0] != s[1])
1387     throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1388     }
1389     else if (getDataPointRank()==4) {
1390     if(!(s[0] == s[2] && s[1] == s[3]))
1391     throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1392     }
1393     else {
1394     throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1395     }
1396     Data ev(0.,getDataPointShape(),getFunctionSpace());
1397     ev.typeMatchRight(*this);
1398     m_data->symmetric(ev.m_data.get());
1399     return ev;
1400     }
1401    
1402     Data
1403     Data::nonsymmetric() const
1404     {
1405     #if defined DOPROF
1406     profData->unary++;
1407     #endif
1408     // check input
1409     DataArrayView::ShapeType s=getDataPointShape();
1410     if (getDataPointRank()==2) {
1411     if(s[0] != s[1])
1412     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1413     DataArrayView::ShapeType ev_shape;
1414     ev_shape.push_back(s[0]);
1415     ev_shape.push_back(s[1]);
1416     Data ev(0.,ev_shape,getFunctionSpace());
1417     ev.typeMatchRight(*this);
1418     m_data->nonsymmetric(ev.m_data.get());
1419     return ev;
1420     }
1421     else if (getDataPointRank()==4) {
1422     if(!(s[0] == s[2] && s[1] == s[3]))
1423     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1424     DataArrayView::ShapeType ev_shape;
1425     ev_shape.push_back(s[0]);
1426     ev_shape.push_back(s[1]);
1427     ev_shape.push_back(s[2]);
1428     ev_shape.push_back(s[3]);
1429     Data ev(0.,ev_shape,getFunctionSpace());
1430     ev.typeMatchRight(*this);
1431     m_data->nonsymmetric(ev.m_data.get());
1432     return ev;
1433     }
1434     else {
1435     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1436     }
1437     }
1438    
1439     Data
1440 gross 800 Data::trace(int axis_offset) const
1441 ksteube 775 {
1442     #if defined DOPROF
1443     profData->unary++;
1444     #endif
1445     DataArrayView::ShapeType s=getDataPointShape();
1446     if (getDataPointRank()==2) {
1447     DataArrayView::ShapeType ev_shape;
1448     Data ev(0.,ev_shape,getFunctionSpace());
1449     ev.typeMatchRight(*this);
1450 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1451 ksteube 775 return ev;
1452     }
1453     if (getDataPointRank()==3) {
1454     DataArrayView::ShapeType ev_shape;
1455     if (axis_offset==0) {
1456     int s2=s[2];
1457     ev_shape.push_back(s2);
1458     }
1459     else if (axis_offset==1) {
1460     int s0=s[0];
1461     ev_shape.push_back(s0);
1462     }
1463     Data ev(0.,ev_shape,getFunctionSpace());
1464     ev.typeMatchRight(*this);
1465 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1466 ksteube 775 return ev;
1467     }
1468     if (getDataPointRank()==4) {
1469     DataArrayView::ShapeType ev_shape;
1470     if (axis_offset==0) {
1471     ev_shape.push_back(s[2]);
1472     ev_shape.push_back(s[3]);
1473     }
1474     else if (axis_offset==1) {
1475     ev_shape.push_back(s[0]);
1476     ev_shape.push_back(s[3]);
1477     }
1478     else if (axis_offset==2) {
1479     ev_shape.push_back(s[0]);
1480     ev_shape.push_back(s[1]);
1481     }
1482     Data ev(0.,ev_shape,getFunctionSpace());
1483     ev.typeMatchRight(*this);
1484 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1485 ksteube 775 return ev;
1486     }
1487     else {
1488 gross 800 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1489 ksteube 775 }
1490     }
1491    
1492     Data
1493     Data::transpose(int axis_offset) const
1494     {
1495 gross 800 #if defined DOPROF
1496     profData->unary++;
1497     #endif
1498 ksteube 775 DataArrayView::ShapeType s=getDataPointShape();
1499     DataArrayView::ShapeType ev_shape;
1500     // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1501     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1502     int rank=getDataPointRank();
1503     if (axis_offset<0 || axis_offset>rank) {
1504     throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1505     }
1506     for (int i=0; i<rank; i++) {
1507     int index = (axis_offset+i)%rank;
1508     ev_shape.push_back(s[index]); // Append to new shape
1509     }
1510     Data ev(0.,ev_shape,getFunctionSpace());
1511     ev.typeMatchRight(*this);
1512     m_data->transpose(ev.m_data.get(), axis_offset);
1513     return ev;
1514 jgs 123 }
1515    
1516 gross 576 Data
1517     Data::eigenvalues() const
1518     {
1519     #if defined DOPROF
1520     profData->unary++;
1521     #endif
1522     // check input
1523     DataArrayView::ShapeType s=getDataPointShape();
1524     if (getDataPointRank()!=2)
1525     throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1526     if(s[0] != s[1])
1527     throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1528     // create return
1529     DataArrayView::ShapeType ev_shape(1,s[0]);
1530     Data ev(0.,ev_shape,getFunctionSpace());
1531     ev.typeMatchRight(*this);
1532     m_data->eigenvalues(ev.m_data.get());
1533     return ev;
1534     }
1535    
1536 jgs 121 const boost::python::tuple
1537 gross 576 Data::eigenvalues_and_eigenvectors(const double tol) const
1538     {
1539     #if defined DOPROF
1540     profData->unary++;
1541     #endif
1542     DataArrayView::ShapeType s=getDataPointShape();
1543     if (getDataPointRank()!=2)
1544     throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1545     if(s[0] != s[1])
1546     throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1547     // create return
1548     DataArrayView::ShapeType ev_shape(1,s[0]);
1549     Data ev(0.,ev_shape,getFunctionSpace());
1550     ev.typeMatchRight(*this);
1551     DataArrayView::ShapeType V_shape(2,s[0]);
1552     Data V(0.,V_shape,getFunctionSpace());
1553     V.typeMatchRight(*this);
1554     m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1555     return make_tuple(boost::python::object(ev),boost::python::object(V));
1556     }
1557    
1558     const boost::python::tuple
1559 jgs 121 Data::mindp() const
1560     {
1561 jgs 148 // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1562     // abort (for unknown reasons) if there are openmp directives with it in the
1563     // surrounding function
1564    
1565     int SampleNo;
1566     int DataPointNo;
1567 bcumming 782 int ProcNo;
1568     calc_mindp(ProcNo,SampleNo,DataPointNo);
1569     return make_tuple(ProcNo,SampleNo,DataPointNo);
1570 jgs 148 }
1571    
1572     void
1573 bcumming 782 Data::calc_mindp( int& ProcNo,
1574     int& SampleNo,
1575     int& DataPointNo) const
1576 jgs 148 {
1577     int i,j;
1578     int lowi=0,lowj=0;
1579     double min=numeric_limits<double>::max();
1580    
1581 jgs 121 Data temp=minval();
1582    
1583     int numSamples=temp.getNumSamples();
1584     int numDPPSample=temp.getNumDataPointsPerSample();
1585    
1586 jgs 148 double next,local_min;
1587     int local_lowi,local_lowj;
1588 jgs 121
1589 jgs 148 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1590     {
1591     local_min=min;
1592     #pragma omp for private(i,j) schedule(static)
1593     for (i=0; i<numSamples; i++) {
1594     for (j=0; j<numDPPSample; j++) {
1595     next=temp.getDataPoint(i,j)();
1596     if (next<local_min) {
1597     local_min=next;
1598     local_lowi=i;
1599     local_lowj=j;
1600     }
1601 jgs 121 }
1602     }
1603 jgs 148 #pragma omp critical
1604     if (local_min<min) {
1605     min=local_min;
1606     lowi=local_lowi;
1607     lowj=local_lowj;
1608     }
1609 jgs 121 }
1610    
1611 bcumming 782 #ifdef PASO_MPI
1612     // determine the processor on which the minimum occurs
1613     next = temp.getDataPoint(lowi,lowj)();
1614     int lowProc = 0;
1615     double *globalMins = new double[get_MPISize()+1];
1616     int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1617    
1618     if( get_MPIRank()==0 ){
1619     next = globalMins[lowProc];
1620     for( i=1; i<get_MPISize(); i++ )
1621     if( next>globalMins[i] ){
1622     lowProc = i;
1623     next = globalMins[i];
1624     }
1625     }
1626     MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1627    
1628     delete [] globalMins;
1629     ProcNo = lowProc;
1630 bcumming 790 #else
1631     ProcNo = 0;
1632 bcumming 782 #endif
1633 jgs 148 SampleNo = lowi;
1634     DataPointNo = lowj;
1635 jgs 121 }
1636    
1637 jgs 104 void
1638     Data::saveDX(std::string fileName) const
1639     {
1640 jgs 153 boost::python::dict args;
1641     args["data"]=boost::python::object(this);
1642     getDomain().saveDX(fileName,args);
1643 jgs 104 return;
1644     }
1645    
1646 jgs 110 void
1647     Data::saveVTK(std::string fileName) const
1648     {
1649 jgs 153 boost::python::dict args;
1650     args["data"]=boost::python::object(this);
1651     getDomain().saveVTK(fileName,args);
1652 jgs 110 return;
1653     }
1654    
1655 jgs 102 Data&
1656     Data::operator+=(const Data& right)
1657     {
1658 gross 783 if (isProtected()) {
1659     throw DataException("Error - attempt to update protected Data object.");
1660     }
1661 jgs 94 binaryOp(right,plus<double>());
1662     return (*this);
1663     }
1664    
1665 jgs 102 Data&
1666     Data::operator+=(const boost::python::object& right)
1667 jgs 94 {
1668 gross 854 Data tmp(right,getFunctionSpace(),false);
1669     binaryOp(tmp,plus<double>());
1670 jgs 94 return (*this);
1671     }
1672    
1673 jgs 102 Data&
1674     Data::operator-=(const Data& right)
1675 jgs 94 {
1676 gross 783 if (isProtected()) {
1677     throw DataException("Error - attempt to update protected Data object.");
1678     }
1679 jgs 94 binaryOp(right,minus<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,minus<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,multiplies<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,multiplies<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,divides<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,divides<double>());
1724 jgs 94 return (*this);
1725     }
1726    
1727 jgs 102 Data
1728 gross 699 Data::rpowO(const boost::python::object& left) const
1729     {
1730     Data left_d(left,*this);
1731     return left_d.powD(*this);
1732     }
1733    
1734     Data
1735 jgs 102 Data::powO(const boost::python::object& right) const
1736 jgs 94 {
1737 gross 854 Data tmp(right,getFunctionSpace(),false);
1738     return powD(tmp);
1739 jgs 94 }
1740    
1741 jgs 102 Data
1742     Data::powD(const Data& right) const
1743 jgs 94 {
1744     Data result;
1745 gross 854 if (getDataPointRank()<right.getDataPointRank()) {
1746     result.copy(right);
1747     result.binaryOp(*this,escript::rpow);
1748     } else {
1749     result.copy(*this);
1750     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1751     }
1752 jgs 94 return result;
1753     }
1754    
1755 bcumming 751
1756 jgs 94 //
1757 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1758 jgs 102 Data
1759     escript::operator+(const Data& left, const Data& right)
1760 jgs 94 {
1761     Data result;
1762     //
1763     // perform a deep copy
1764 gross 854 if (left.getDataPointRank()<right.getDataPointRank()) {
1765     result.copy(right);
1766     result+=left;
1767     } else {
1768     result.copy(left);
1769     result+=right;
1770     }
1771 jgs 94 return result;
1772     }
1773    
1774     //
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=right.neg();
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.copy(right);
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=right.oneOver();
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 boost::python::object& right)
1832 jgs 94 {
1833 gross 854 return left+Data(right,left.getFunctionSpace(),false);
1834 jgs 94 }
1835    
1836     //
1837 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1838 jgs 102 Data
1839     escript::operator-(const Data& left, const boost::python::object& right)
1840 jgs 94 {
1841 gross 854 return left-Data(right,left.getFunctionSpace(),false);
1842 jgs 94 }
1843    
1844     //
1845 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1846 jgs 102 Data
1847     escript::operator*(const Data& left, const boost::python::object& right)
1848 jgs 94 {
1849 gross 854 return left*Data(right,left.getFunctionSpace(),false);
1850 jgs 94 }
1851    
1852     //
1853 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1854 jgs 102 Data
1855     escript::operator/(const Data& left, const boost::python::object& right)
1856 jgs 94 {
1857 gross 854 return left/Data(right,left.getFunctionSpace(),false);
1858 jgs 94 }
1859    
1860     //
1861 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1862 jgs 102 Data
1863     escript::operator+(const boost::python::object& left, const Data& right)
1864 jgs 94 {
1865 gross 854 return Data(left,right.getFunctionSpace(),false)+right;
1866 jgs 94 }
1867    
1868     //
1869 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1870 jgs 102 Data
1871     escript::operator-(const boost::python::object& left, const Data& right)
1872 jgs 94 {
1873 gross 854 return Data(left,right.getFunctionSpace(),false)-right;
1874 jgs 94 }
1875    
1876     //
1877 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1878 jgs 102 Data
1879     escript::operator*(const boost::python::object& left, const Data& right)
1880 jgs 94 {
1881 gross 854 return Data(left,right.getFunctionSpace(),false)*right;
1882 jgs 94 }
1883    
1884     //
1885 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1886 jgs 102 Data
1887     escript::operator/(const boost::python::object& left, const Data& right)
1888 jgs 94 {
1889 gross 854 return Data(left,right.getFunctionSpace(),false)/right;
1890 jgs 94 }
1891    
1892     //
1893 jgs 102 //bool escript::operator==(const Data& left, const Data& right)
1894     //{
1895     // /*
1896     // NB: this operator does very little at this point, and isn't to
1897     // be relied on. Requires further implementation.
1898     // */
1899     //
1900     // bool ret;
1901     //
1902     // if (left.isEmpty()) {
1903     // if(!right.isEmpty()) {
1904     // ret = false;
1905     // } else {
1906     // ret = true;
1907     // }
1908     // }
1909     //
1910     // if (left.isConstant()) {
1911     // if(!right.isConstant()) {
1912     // ret = false;
1913     // } else {
1914     // ret = true;
1915     // }
1916     // }
1917     //
1918     // if (left.isTagged()) {
1919     // if(!right.isTagged()) {
1920     // ret = false;
1921     // } else {
1922     // ret = true;
1923     // }
1924     // }
1925     //
1926     // if (left.isExpanded()) {
1927     // if(!right.isExpanded()) {
1928     // ret = false;
1929     // } else {
1930     // ret = true;
1931     // }
1932     // }
1933     //
1934     // return ret;
1935     //}
1936    
1937 bcumming 751 /* TODO */
1938     /* global reduction */
1939 jgs 102 Data
1940     Data::getItem(const boost::python::object& key) const
1941 jgs 94 {
1942 jgs 102 const DataArrayView& view=getPointDataView();
1943 jgs 94
1944 jgs 102 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1945 jgs 94
1946 jgs 102 if (slice_region.size()!=view.getRank()) {
1947     throw DataException("Error - slice size does not match Data rank.");
1948 jgs 94 }
1949    
1950 jgs 102 return getSlice(slice_region);
1951 jgs 94 }
1952    
1953 bcumming 751 /* TODO */
1954     /* global reduction */
1955 jgs 94 Data
1956 jgs 102 Data::getSlice(const DataArrayView::RegionType& region) const
1957 jgs 94 {
1958 jgs 151 #if defined DOPROF
1959 jgs 123 profData->slicing++;
1960 jgs 151 #endif
1961 jgs 102 return Data(*this,region);
1962 jgs 94 }
1963    
1964 bcumming 751 /* TODO */
1965     /* global reduction */
1966 jgs 94 void
1967 jgs 102 Data::setItemO(const boost::python::object& key,
1968     const boost::python::object& value)
1969 jgs 94 {
1970 jgs 102 Data tempData(value,getFunctionSpace());
1971     setItemD(key,tempData);
1972     }
1973    
1974 bcumming 751 /* TODO */
1975     /* global reduction */
1976 jgs 102 void
1977     Data::setItemD(const boost::python::object& key,
1978     const Data& value)
1979     {
1980 jgs 94 const DataArrayView& view=getPointDataView();
1981 jgs 123
1982 jgs 94 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1983     if (slice_region.size()!=view.getRank()) {
1984     throw DataException("Error - slice size does not match Data rank.");
1985     }
1986 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
1987     setSlice(Data(value,getFunctionSpace()),slice_region);
1988     } else {
1989     setSlice(value,slice_region);
1990     }
1991 jgs 94 }
1992    
1993 bcumming 751 /* TODO */
1994     /* global reduction */
1995 jgs 94 void
1996 jgs 102 Data::setSlice(const Data& value,
1997     const DataArrayView::RegionType& region)
1998 jgs 94 {
1999 gross 783 if (isProtected()) {
2000     throw DataException("Error - attempt to update protected Data object.");
2001     }
2002 jgs 151 #if defined DOPROF
2003 jgs 123 profData->slicing++;
2004 jgs 151 #endif
2005 jgs 102 Data tempValue(value);
2006     typeMatchLeft(tempValue);
2007     typeMatchRight(tempValue);
2008     m_data->setSlice(tempValue.m_data.get(),region);
2009     }
2010    
2011     void
2012     Data::typeMatchLeft(Data& right) const
2013     {
2014     if (isExpanded()){
2015     right.expand();
2016     } else if (isTagged()) {
2017     if (right.isConstant()) {
2018     right.tag();
2019     }
2020     }
2021     }
2022    
2023     void
2024     Data::typeMatchRight(const Data& right)
2025     {
2026 jgs 94 if (isTagged()) {
2027     if (right.isExpanded()) {
2028     expand();
2029     }
2030     } else if (isConstant()) {
2031     if (right.isExpanded()) {
2032     expand();
2033     } else if (right.isTagged()) {
2034     tag();
2035     }
2036     }
2037     }
2038    
2039 bcumming 751 /* TODO */
2040     /* global reduction */
2041 jgs 94 void
2042     Data::setTaggedValue(int tagKey,
2043     const boost::python::object& value)
2044     {
2045 gross 783 if (isProtected()) {
2046     throw DataException("Error - attempt to update protected Data object.");
2047     }
2048 jgs 94 //
2049     // Ensure underlying data object is of type DataTagged
2050     tag();
2051    
2052     if (!isTagged()) {
2053     throw DataException("Error - DataTagged conversion failed!!");
2054     }
2055    
2056     //
2057     // Construct DataArray from boost::python::object input value
2058     DataArray valueDataArray(value);
2059    
2060     //
2061     // Call DataAbstract::setTaggedValue
2062     m_data->setTaggedValue(tagKey,valueDataArray.getView());
2063     }
2064    
2065 bcumming 751 /* TODO */
2066     /* global reduction */
2067 jgs 110 void
2068 jgs 121 Data::setTaggedValueFromCPP(int tagKey,
2069     const DataArrayView& value)
2070     {
2071 gross 783 if (isProtected()) {
2072     throw DataException("Error - attempt to update protected Data object.");
2073     }
2074 jgs 121 //
2075     // Ensure underlying data object is of type DataTagged
2076     tag();
2077    
2078     if (!isTagged()) {
2079     throw DataException("Error - DataTagged conversion failed!!");
2080     }
2081    
2082     //
2083     // Call DataAbstract::setTaggedValue
2084     m_data->setTaggedValue(tagKey,value);
2085     }
2086    
2087 bcumming 751 /* TODO */
2088     /* global reduction */
2089 jgs 149 int
2090     Data::getTagNumber(int dpno)
2091     {
2092     return m_data->getTagNumber(dpno);
2093     }
2094    
2095 bcumming 751 /* TODO */
2096     /* global reduction */
2097 jgs 121 void
2098 jgs 110 Data::setRefValue(int ref,
2099     const boost::python::numeric::array& value)
2100     {
2101 gross 783 if (isProtected()) {
2102     throw DataException("Error - attempt to update protected Data object.");
2103     }
2104 jgs 110 //
2105     // Construct DataArray from boost::python::object input value
2106     DataArray valueDataArray(value);
2107    
2108     //
2109     // Call DataAbstract::setRefValue
2110     m_data->setRefValue(ref,valueDataArray);
2111     }
2112    
2113 bcumming 751 /* TODO */
2114     /* global reduction */
2115 jgs 110 void
2116     Data::getRefValue(int ref,
2117     boost::python::numeric::array& value)
2118     {
2119     //
2120     // Construct DataArray for boost::python::object return value
2121     DataArray valueDataArray(value);
2122    
2123     //
2124     // Load DataArray with values from data-points specified by ref
2125     m_data->getRefValue(ref,valueDataArray);
2126    
2127     //
2128     // Load values from valueDataArray into return numarray
2129    
2130     // extract the shape of the numarray
2131     int rank = value.getrank();
2132     DataArrayView::ShapeType shape;
2133     for (int i=0; i < rank; i++) {
2134     shape.push_back(extract<int>(value.getshape()[i]));
2135     }
2136    
2137     // and load the numarray with the data from the DataArray
2138     DataArrayView valueView = valueDataArray.getView();
2139    
2140     if (rank==0) {
2141 jgs 126 boost::python::numeric::array temp_numArray(valueView());
2142     value = temp_numArray;
2143 jgs 110 }
2144     if (rank==1) {
2145     for (int i=0; i < shape[0]; i++) {
2146     value[i] = valueView(i);
2147     }
2148     }
2149     if (rank==2) {
2150 jgs 126 for (int i=0; i < shape[0]; i++) {
2151     for (int j=0; j < shape[1]; j++) {
2152     value[i][j] = valueView(i,j);
2153     }
2154     }
2155 jgs 110 }
2156     if (rank==3) {
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     value[i][j][k] = valueView(i,j,k);
2161     }
2162     }
2163     }
2164 jgs 110 }
2165     if (rank==4) {
2166 jgs 126 for (int i=0; i < shape[0]; i++) {
2167     for (int j=0; j < shape[1]; j++) {
2168     for (int k=0; k < shape[2]; k++) {
2169     for (int l=0; l < shape[3]; l++) {
2170     value[i][j][k][l] = valueView(i,j,k,l);
2171     }
2172     }
2173     }
2174     }
2175 jgs 110 }
2176    
2177     }
2178    
2179 jgs 94 void
2180 jgs 119 Data::archiveData(const std::string fileName)
2181     {
2182     cout << "Archiving Data object to: " << fileName << endl;
2183    
2184     //
2185     // Determine type of this Data object
2186     int dataType = -1;
2187    
2188     if (isEmpty()) {
2189     dataType = 0;
2190     cout << "\tdataType: DataEmpty" << endl;
2191     }
2192     if (isConstant()) {
2193     dataType = 1;
2194     cout << "\tdataType: DataConstant" << endl;
2195     }
2196     if (isTagged()) {
2197     dataType = 2;
2198     cout << "\tdataType: DataTagged" << endl;
2199     }
2200     if (isExpanded()) {
2201     dataType = 3;
2202     cout << "\tdataType: DataExpanded" << endl;
2203     }
2204 jgs 123
2205 jgs 119 if (dataType == -1) {
2206     throw DataException("archiveData Error: undefined dataType");
2207     }
2208    
2209     //
2210     // Collect data items common to all Data types
2211     int noSamples = getNumSamples();
2212     int noDPPSample = getNumDataPointsPerSample();
2213     int functionSpaceType = getFunctionSpace().getTypeCode();
2214     int dataPointRank = getDataPointRank();
2215     int dataPointSize = getDataPointSize();
2216     int dataLength = getLength();
2217     DataArrayView::ShapeType dataPointShape = getDataPointShape();
2218 woo409 757 vector<int> referenceNumbers(noSamples);
2219 jgs 119 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2220     referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
2221     }
2222 woo409 757 vector<int> tagNumbers(noSamples);
2223 jgs 119 if (isTagged()) {
2224     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2225     tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
2226     }
2227     }
2228    
2229     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2230     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2231     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2232    
2233     //
2234     // Flatten Shape to an array of integers suitable for writing to file
2235     int flatShape[4] = {0,0,0,0};
2236     cout << "\tshape: < ";
2237     for (int dim=0; dim<dataPointRank; dim++) {
2238     flatShape[dim] = dataPointShape[dim];
2239     cout << dataPointShape[dim] << " ";
2240     }
2241     cout << ">" << endl;
2242    
2243     //
2244 jgs 123 // Open archive file
2245 jgs 119 ofstream archiveFile;
2246     archiveFile.open(fileName.data(), ios::out);
2247    
2248     if (!archiveFile.good()) {
2249     throw DataException("archiveData Error: problem opening archive file");
2250     }
2251    
2252 jgs 123 //
2253     // Write common data items to archive file
2254 jgs 119 archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
2255     archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
2256     archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2257     archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2258     archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2259     archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2260     archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
2261     for (int dim = 0; dim < 4; dim++) {
2262     archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2263     }
2264     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2265     archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2266     }
2267     if (isTagged()) {
2268     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2269     archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2270     }
2271     }
2272    
2273     if (!archiveFile.good()) {
2274     throw DataException("archiveData Error: problem writing to archive file");
2275     }
2276    
2277     //
2278 jgs 123 // Archive underlying data values for each Data type
2279     int noValues;
2280 jgs 119 switch (dataType) {
2281     case 0:
2282     // DataEmpty
2283 jgs 123 noValues = 0;
2284     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2285     cout << "\tnoValues: " << noValues << endl;
2286 jgs 119 break;
2287     case 1:
2288     // DataConstant
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 2:
2297     // DataTagged
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     case 3:
2306     // DataExpanded
2307 jgs 123 noValues = m_data->getLength();
2308     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2309     cout << "\tnoValues: " << noValues << endl;
2310     if (m_data->archiveData(archiveFile,noValues)) {
2311     throw DataException("archiveData Error: problem writing data to archive file");
2312     }
2313 jgs 119 break;
2314     }
2315    
2316 jgs 123 if (!archiveFile.good()) {
2317     throw DataException("archiveData Error: problem writing data to archive file");
2318     }
2319    
2320     //
2321     // Close archive file
2322     archiveFile.close();
2323    
2324     if (!archiveFile.good()) {
2325     throw DataException("archiveData Error: problem closing archive file");
2326     }
2327    
2328 jgs 119 }
2329    
2330     void
2331     Data::extractData(const std::string fileName,
2332     const FunctionSpace& fspace)
2333     {
2334     //
2335     // Can only extract Data to an object which is initially DataEmpty
2336     if (!isEmpty()) {
2337     throw DataException("extractData Error: can only extract to DataEmpty object");
2338     }
2339    
2340     cout << "Extracting Data object from: " << fileName << endl;
2341    
2342     int dataType;
2343     int noSamples;
2344     int noDPPSample;
2345     int functionSpaceType;
2346     int dataPointRank;
2347     int dataPointSize;
2348     int dataLength;
2349     DataArrayView::ShapeType dataPointShape;
2350     int flatShape[4];
2351    
2352     //
2353 jgs 123 // Open the archive file
2354 jgs 119 ifstream archiveFile;
2355     archiveFile.open(fileName.data(), ios::in);
2356    
2357     if (!archiveFile.good()) {
2358     throw DataException("extractData Error: problem opening archive file");
2359     }
2360    
2361 jgs 123 //
2362     // Read common data items from archive file
2363 jgs 119 archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
2364     archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2365     archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2366     archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2367     archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2368     archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2369     archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
2370     for (int dim = 0; dim < 4; dim++) {
2371     archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2372     if (flatShape[dim]>0) {
2373     dataPointShape.push_back(flatShape[dim]);
2374     }
2375     }
2376 woo409 757 vector<int> referenceNumbers(noSamples);
2377 jgs 119 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2378     archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2379     }
2380 woo409 757 vector<int> tagNumbers(noSamples);
2381 jgs 119 if (dataType==2) {
2382     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2383     archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2384     }
2385     }
2386    
2387     if (!archiveFile.good()) {
2388     throw DataException("extractData Error: problem reading from archive file");
2389     }
2390    
2391 jgs 123 //
2392     // Verify the values just read from the archive file
2393 jgs 119 switch (dataType) {
2394     case 0:
2395     cout << "\tdataType: DataEmpty" << endl;
2396     break;
2397     case 1:
2398     cout << "\tdataType: DataConstant" << endl;
2399     break;
2400     case 2:
2401     cout << "\tdataType: DataTagged" << endl;
2402     break;
2403     case 3:
2404     cout << "\tdataType: DataExpanded" << endl;
2405     break;
2406     default:
2407     throw DataException("extractData Error: undefined dataType read from archive file");
2408     break;
2409     }
2410    
2411     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2412     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2413     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2414     cout << "\tshape: < ";
2415     for (int dim = 0; dim < dataPointRank; dim++) {
2416     cout << dataPointShape[dim] << " ";
2417     }
2418     cout << ">" << endl;
2419    
2420     //
2421     // Verify that supplied FunctionSpace object is compatible with this Data object.
2422     if ( (fspace.getTypeCode()!=functionSpaceType) ||
2423     (fspace.getNumSamples()!=noSamples) ||
2424     (fspace.getNumDPPSample()!=noDPPSample)
2425     ) {
2426     throw DataException("extractData Error: incompatible FunctionSpace");
2427     }
2428     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2429     if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2430     throw DataException("extractData Error: incompatible FunctionSpace");
2431     }
2432     }
2433     if (dataType==2) {
2434     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2435     if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2436     throw DataException("extractData Error: incompatible FunctionSpace");
2437     }
2438     }
2439     }
2440    
2441     //
2442     // Construct a DataVector to hold underlying data values
2443     DataVector dataVec(dataLength);
2444    
2445     //
2446     // Load this DataVector with the appropriate values
2447 jgs 123 int noValues;
2448     archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2449     cout << "\tnoValues: " << noValues << endl;
2450 jgs 119 switch (dataType) {
2451     case 0:
2452     // DataEmpty
2453 jgs 123 if (noValues != 0) {
2454     throw DataException("extractData Error: problem reading data from archive file");
2455     }
2456 jgs 119 break;
2457     case 1:
2458     // DataConstant
2459 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2460     throw DataException("extractData Error: problem reading data from archive file");
2461     }
2462 jgs 119 break;
2463     case 2:
2464     // DataTagged
2465 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2466     throw DataException("extractData Error: problem reading data from archive file");
2467     }
2468 jgs 119 break;
2469     case 3:
2470     // DataExpanded
2471 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2472     throw DataException("extractData Error: problem reading data from archive file");
2473     }
2474 jgs 119 break;
2475     }
2476    
2477 jgs 123 if (!archiveFile.good()) {
2478     throw DataException("extractData Error: problem reading from archive file");
2479     }
2480    
2481 jgs 119 //
2482 jgs 123 // Close archive file
2483     archiveFile.close();
2484    
2485     if (!archiveFile.good()) {
2486     throw DataException("extractData Error: problem closing archive file");
2487     }
2488    
2489     //
2490 jgs 119 // Construct an appropriate Data object
2491     DataAbstract* tempData;
2492     switch (dataType) {
2493     case 0:
2494     // DataEmpty
2495     tempData=new DataEmpty();