/[escript]/branches/intelc_win32/escript/src/Data.cpp
ViewVC logotype

Annotation of /branches/intelc_win32/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26