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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 757 - (hide annotations)
Mon Jun 26 13:12:56 2006 UTC (13 years, 2 months ago) by woo409
File size: 54977 byte(s)
+ Merge of intelc_win32 branch (revision 741:755) with trunk. Tested on iVEC altix (run_tests and py_tests all pass)

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26