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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 711 - (hide annotations)
Wed Apr 26 22:39:51 2006 UTC (13 years, 8 months ago) by gross
File size: 53373 byte(s)
tests pass now on gcc: some tests did not take round-off errors into consideration
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 94
15 jgs 474 #include "Data.h"
16 jgs 94
17 jgs 480 #include "DataExpanded.h"
18     #include "DataConstant.h"
19     #include "DataTagged.h"
20     #include "DataEmpty.h"
21     #include "DataArray.h"
22     #include "DataArrayView.h"
23     #include "DataProf.h"
24     #include "FunctionSpaceFactory.h"
25     #include "AbstractContinuousDomain.h"
26     #include "UnaryFuncs.h"
27    
28 jgs 119 #include <fstream>
29 jgs 94 #include <algorithm>
30     #include <vector>
31     #include <functional>
32     #include <math.h>
33    
34 jgs 153 #include <boost/python/dict.hpp>
35 jgs 94 #include <boost/python/extract.hpp>
36     #include <boost/python/long.hpp>
37    
38     using namespace std;
39     using namespace boost::python;
40     using namespace boost;
41     using namespace escript;
42    
43 jgs 151 #if defined DOPROF
44 jgs 123 //
45     // global table of profiling data for all Data objects
46     DataProf dataProfTable;
47 jgs 151 #endif
48 jgs 123
49 jgs 94 Data::Data()
50     {
51     //
52     // Default data is type DataEmpty
53     DataAbstract* temp=new DataEmpty();
54 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
55     m_data=temp_data;
56 jgs 151 #if defined DOPROF
57 jgs 123 // create entry in global profiling table for this object
58     profData = dataProfTable.newData();
59 jgs 151 #endif
60 jgs 94 }
61    
62     Data::Data(double value,
63     const tuple& shape,
64     const FunctionSpace& what,
65     bool expanded)
66     {
67     DataArrayView::ShapeType dataPointShape;
68     for (int i = 0; i < shape.attr("__len__")(); ++i) {
69     dataPointShape.push_back(extract<const int>(shape[i]));
70     }
71     DataArray temp(dataPointShape,value);
72     initialise(temp.getView(),what,expanded);
73 jgs 151 #if defined DOPROF
74 jgs 123 // create entry in global profiling table for this object
75     profData = dataProfTable.newData();
76 jgs 151 #endif
77 jgs 94 }
78    
79     Data::Data(double value,
80     const DataArrayView::ShapeType& dataPointShape,
81     const FunctionSpace& what,
82     bool expanded)
83     {
84     DataArray temp(dataPointShape,value);
85     pair<int,int> dataShape=what.getDataShape();
86     initialise(temp.getView(),what,expanded);
87 jgs 151 #if defined DOPROF
88 jgs 123 // create entry in global profiling table for this object
89     profData = dataProfTable.newData();
90 jgs 151 #endif
91 jgs 94 }
92    
93 jgs 102 Data::Data(const Data& inData)
94 jgs 94 {
95 jgs 102 m_data=inData.m_data;
96 jgs 151 #if defined DOPROF
97 jgs 123 // create entry in global profiling table for this object
98     profData = dataProfTable.newData();
99 jgs 151 #endif
100 jgs 94 }
101    
102     Data::Data(const Data& inData,
103     const DataArrayView::RegionType& region)
104     {
105     //
106 jgs 102 // Create Data which is a slice of another Data
107     DataAbstract* tmp = inData.m_data->getSlice(region);
108     shared_ptr<DataAbstract> temp_data(tmp);
109     m_data=temp_data;
110 jgs 151 #if defined DOPROF
111 jgs 123 // create entry in global profiling table for this object
112     profData = dataProfTable.newData();
113 jgs 151 #endif
114 jgs 94 }
115    
116     Data::Data(const Data& inData,
117     const FunctionSpace& functionspace)
118     {
119 gross 711 #if defined DOPROF
120     // create entry in global profiling table for this object
121     profData = dataProfTable.newData();
122     #endif
123 jgs 94 if (inData.getFunctionSpace()==functionspace) {
124     m_data=inData.m_data;
125     } else {
126 gross 711 #if defined DOPROF
127     profData->interpolate++;
128     #endif
129 jgs 94 Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
130 jgs 123 // Note: Must use a reference or pointer to a derived object
131 jgs 94 // in order to get polymorphic behaviour. Shouldn't really
132     // be able to create an instance of AbstractDomain but that was done
133 jgs 123 // as a boost:python work around which may no longer be required.
134 jgs 94 const AbstractDomain& inDataDomain=inData.getDomain();
135     if (inDataDomain==functionspace.getDomain()) {
136     inDataDomain.interpolateOnDomain(tmp,inData);
137     } else {
138     inDataDomain.interpolateACross(tmp,inData);
139     }
140     m_data=tmp.m_data;
141     }
142     }
143    
144     Data::Data(const DataTagged::TagListType& tagKeys,
145     const DataTagged::ValueListType & values,
146     const DataArrayView& defaultValue,
147     const FunctionSpace& what,
148     bool expanded)
149     {
150     DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
151 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
152     m_data=temp_data;
153 jgs 94 if (expanded) {
154     expand();
155     }
156 jgs 151 #if defined DOPROF
157 jgs 123 // create entry in global profiling table for this object
158     profData = dataProfTable.newData();
159 jgs 151 #endif
160 jgs 94 }
161    
162     Data::Data(const numeric::array& value,
163     const FunctionSpace& what,
164     bool expanded)
165     {
166     initialise(value,what,expanded);
167 jgs 151 #if defined DOPROF
168 jgs 123 // create entry in global profiling table for this object
169     profData = dataProfTable.newData();
170 jgs 151 #endif
171 jgs 94 }
172    
173     Data::Data(const DataArrayView& value,
174     const FunctionSpace& what,
175     bool expanded)
176     {
177     initialise(value,what,expanded);
178 jgs 151 #if defined DOPROF
179 jgs 123 // create entry in global profiling table for this object
180     profData = dataProfTable.newData();
181 jgs 151 #endif
182 jgs 94 }
183    
184     Data::Data(const object& value,
185     const FunctionSpace& what,
186     bool expanded)
187     {
188     numeric::array asNumArray(value);
189     initialise(asNumArray,what,expanded);
190 jgs 151 #if defined DOPROF
191 jgs 123 // create entry in global profiling table for this object
192     profData = dataProfTable.newData();
193 jgs 151 #endif
194 jgs 94 }
195    
196     Data::Data(const object& value,
197     const Data& other)
198     {
199     //
200     // Create DataConstant using the given value and all other parameters
201     // copied from other. If value is a rank 0 object this Data
202     // will assume the point data shape of other.
203     DataArray temp(value);
204     if (temp.getView().getRank()==0) {
205     //
206     // Create a DataArray with the scalar value for all elements
207     DataArray temp2(other.getPointDataView().getShape(),temp.getView()());
208     initialise(temp2.getView(),other.getFunctionSpace(),false);
209     } else {
210     //
211     // Create a DataConstant with the same sample shape as other
212     initialise(temp.getView(),other.getFunctionSpace(),false);
213     }
214 jgs 151 #if defined DOPROF
215 jgs 123 // create entry in global profiling table for this object
216     profData = dataProfTable.newData();
217 jgs 151 #endif
218 jgs 94 }
219    
220 jgs 151 Data::~Data()
221     {
222    
223     }
224    
225 jgs 94 escriptDataC
226     Data::getDataC()
227     {
228     escriptDataC temp;
229     temp.m_dataPtr=(void*)this;
230     return temp;
231     }
232    
233     escriptDataC
234     Data::getDataC() const
235     {
236     escriptDataC temp;
237     temp.m_dataPtr=(void*)this;
238     return temp;
239     }
240    
241 jgs 121 const boost::python::tuple
242 jgs 94 Data::getShapeTuple() const
243     {
244     const DataArrayView::ShapeType& shape=getDataPointShape();
245     switch(getDataPointRank()) {
246     case 0:
247     return make_tuple();
248     case 1:
249     return make_tuple(long_(shape[0]));
250     case 2:
251     return make_tuple(long_(shape[0]),long_(shape[1]));
252     case 3:
253     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]));
254     case 4:
255     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]),long_(shape[3]));
256     default:
257     throw DataException("Error - illegal Data rank.");
258     }
259     }
260    
261     void
262     Data::copy(const Data& other)
263     {
264     //
265     // Perform a deep copy
266     {
267     DataExpanded* temp=dynamic_cast<DataExpanded*>(other.m_data.get());
268     if (temp!=0) {
269     //
270     // Construct a DataExpanded copy
271     DataAbstract* newData=new DataExpanded(*temp);
272 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
273     m_data=temp_data;
274 jgs 94 return;
275     }
276     }
277     {
278     DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());
279     if (temp!=0) {
280     //
281 jgs 102 // Construct a DataTagged copy
282 jgs 94 DataAbstract* newData=new DataTagged(*temp);
283 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
284     m_data=temp_data;
285 jgs 94 return;
286     }
287     }
288     {
289     DataConstant* temp=dynamic_cast<DataConstant*>(other.m_data.get());
290     if (temp!=0) {
291     //
292     // Construct a DataConstant copy
293     DataAbstract* newData=new DataConstant(*temp);
294 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
295     m_data=temp_data;
296 jgs 94 return;
297     }
298     }
299 jgs 102 {
300     DataEmpty* temp=dynamic_cast<DataEmpty*>(other.m_data.get());
301     if (temp!=0) {
302     //
303     // Construct a DataEmpty copy
304     DataAbstract* newData=new DataEmpty();
305     shared_ptr<DataAbstract> temp_data(newData);
306     m_data=temp_data;
307     return;
308     }
309     }
310 jgs 94 throw DataException("Error - Copy not implemented for this Data type.");
311     }
312    
313     void
314     Data::copyWithMask(const Data& other,
315     const Data& mask)
316     {
317     Data mask1;
318     Data mask2;
319    
320     mask1 = mask.wherePositive();
321     mask2.copy(mask1);
322    
323     mask1 *= other;
324     mask2 *= *this;
325     mask2 = *this - mask2;
326    
327     *this = mask1 + mask2;
328     }
329    
330     bool
331     Data::isExpanded() const
332     {
333     DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
334     return (temp!=0);
335     }
336    
337     bool
338     Data::isTagged() const
339     {
340     DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
341     return (temp!=0);
342     }
343    
344     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 jgs 151 #if defined DOPROF
1107 jgs 123 profData->reduction1++;
1108 jgs 151 #endif
1109 jgs 106 //
1110     // set the initial absolute maximum value to zero
1111 jgs 147 AbsMax abs_max_func;
1112     return algorithm(abs_max_func,0);
1113 jgs 102 }
1114    
1115 jgs 106 double
1116 jgs 117 Data::Linf() const
1117     {
1118 jgs 151 #if defined DOPROF
1119 jgs 123 profData->reduction1++;
1120 jgs 151 #endif
1121 jgs 117 //
1122     // set the initial absolute minimum value to max double
1123 jgs 147 AbsMin abs_min_func;
1124     return algorithm(abs_min_func,numeric_limits<double>::max());
1125 jgs 117 }
1126    
1127     double
1128 jgs 106 Data::sup() const
1129 jgs 102 {
1130 jgs 151 #if defined DOPROF
1131 jgs 123 profData->reduction1++;
1132 jgs 151 #endif
1133 jgs 106 //
1134     // set the initial maximum value to min possible double
1135 jgs 147 FMax fmax_func;
1136     return algorithm(fmax_func,numeric_limits<double>::max()*-1);
1137 jgs 102 }
1138    
1139 jgs 106 double
1140     Data::inf() const
1141 jgs 102 {
1142 jgs 151 #if defined DOPROF
1143 jgs 123 profData->reduction1++;
1144 jgs 151 #endif
1145 jgs 106 //
1146     // set the initial minimum value to max possible double
1147 jgs 147 FMin fmin_func;
1148     return algorithm(fmin_func,numeric_limits<double>::max());
1149 jgs 102 }
1150    
1151     Data
1152 jgs 106 Data::maxval() const
1153 jgs 102 {
1154 jgs 151 #if defined DOPROF
1155 jgs 123 profData->reduction2++;
1156 jgs 151 #endif
1157 jgs 113 //
1158     // set the initial maximum value to min possible double
1159 jgs 147 FMax fmax_func;
1160     return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1161 jgs 102 }
1162    
1163     Data
1164 jgs 106 Data::minval() const
1165 jgs 102 {
1166 jgs 151 #if defined DOPROF
1167 jgs 123 profData->reduction2++;
1168 jgs 151 #endif
1169 jgs 113 //
1170     // set the initial minimum value to max possible double
1171 jgs 147 FMin fmin_func;
1172     return dp_algorithm(fmin_func,numeric_limits<double>::max());
1173 jgs 102 }
1174    
1175 jgs 123 Data
1176     Data::trace() const
1177     {
1178 jgs 151 #if defined DOPROF
1179 jgs 123 profData->reduction2++;
1180 jgs 151 #endif
1181 jgs 147 Trace trace_func;
1182     return dp_algorithm(trace_func,0);
1183 jgs 123 }
1184    
1185     Data
1186     Data::transpose(int axis) const
1187     {
1188 jgs 151 #if defined DOPROF
1189 jgs 123 profData->reduction2++;
1190 jgs 151 #endif
1191 gross 576
1192 jgs 123 // not implemented
1193     throw DataException("Error - Data::transpose not implemented yet.");
1194     return Data();
1195     }
1196    
1197 gross 576 Data
1198     Data::eigenvalues() const
1199     {
1200     #if defined DOPROF
1201     profData->unary++;
1202     #endif
1203     // check input
1204     DataArrayView::ShapeType s=getDataPointShape();
1205     if (getDataPointRank()!=2)
1206     throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1207     if(s[0] != s[1])
1208     throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1209     // create return
1210     DataArrayView::ShapeType ev_shape(1,s[0]);
1211     Data ev(0.,ev_shape,getFunctionSpace());
1212     ev.typeMatchRight(*this);
1213     m_data->eigenvalues(ev.m_data.get());
1214     return ev;
1215     }
1216    
1217 jgs 121 const boost::python::tuple
1218 gross 576 Data::eigenvalues_and_eigenvectors(const double tol) const
1219     {
1220     #if defined DOPROF
1221     profData->unary++;
1222     #endif
1223     DataArrayView::ShapeType s=getDataPointShape();
1224     if (getDataPointRank()!=2)
1225     throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1226     if(s[0] != s[1])
1227     throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1228     // create return
1229     DataArrayView::ShapeType ev_shape(1,s[0]);
1230     Data ev(0.,ev_shape,getFunctionSpace());
1231     ev.typeMatchRight(*this);
1232     DataArrayView::ShapeType V_shape(2,s[0]);
1233     Data V(0.,V_shape,getFunctionSpace());
1234     V.typeMatchRight(*this);
1235     m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1236     return make_tuple(boost::python::object(ev),boost::python::object(V));
1237     }
1238    
1239     const boost::python::tuple
1240 jgs 121 Data::mindp() const
1241     {
1242 jgs 148 // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1243     // abort (for unknown reasons) if there are openmp directives with it in the
1244     // surrounding function
1245    
1246     int SampleNo;
1247     int DataPointNo;
1248    
1249     calc_mindp(SampleNo,DataPointNo);
1250    
1251     return make_tuple(SampleNo,DataPointNo);
1252     }
1253    
1254     void
1255     Data::calc_mindp(int& SampleNo,
1256     int& DataPointNo) const
1257     {
1258     int i,j;
1259     int lowi=0,lowj=0;
1260     double min=numeric_limits<double>::max();
1261    
1262 jgs 121 Data temp=minval();
1263    
1264     int numSamples=temp.getNumSamples();
1265     int numDPPSample=temp.getNumDataPointsPerSample();
1266    
1267 jgs 148 double next,local_min;
1268     int local_lowi,local_lowj;
1269 jgs 121
1270 jgs 148 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1271     {
1272     local_min=min;
1273     #pragma omp for private(i,j) schedule(static)
1274     for (i=0; i<numSamples; i++) {
1275     for (j=0; j<numDPPSample; j++) {
1276     next=temp.getDataPoint(i,j)();
1277     if (next<local_min) {
1278     local_min=next;
1279     local_lowi=i;
1280     local_lowj=j;
1281     }
1282 jgs 121 }
1283     }
1284 jgs 148 #pragma omp critical
1285     if (local_min<min) {
1286     min=local_min;
1287     lowi=local_lowi;
1288     lowj=local_lowj;
1289     }
1290 jgs 121 }
1291    
1292 jgs 148 SampleNo = lowi;
1293     DataPointNo = lowj;
1294 jgs 121 }
1295    
1296 jgs 104 void
1297     Data::saveDX(std::string fileName) const
1298     {
1299 jgs 153 boost::python::dict args;
1300     args["data"]=boost::python::object(this);
1301     getDomain().saveDX(fileName,args);
1302 jgs 104 return;
1303     }
1304    
1305 jgs 110 void
1306     Data::saveVTK(std::string fileName) const
1307     {
1308 jgs 153 boost::python::dict args;
1309     args["data"]=boost::python::object(this);
1310     getDomain().saveVTK(fileName,args);
1311 jgs 110 return;
1312     }
1313    
1314 jgs 102 Data&
1315     Data::operator+=(const Data& right)
1316     {
1317 jgs 151 #if defined DOPROF
1318 jgs 123 profData->binary++;
1319 jgs 151 #endif
1320 jgs 94 binaryOp(right,plus<double>());
1321     return (*this);
1322     }
1323    
1324 jgs 102 Data&
1325     Data::operator+=(const boost::python::object& right)
1326 jgs 94 {
1327 jgs 151 #if defined DOPROF
1328 jgs 123 profData->binary++;
1329 jgs 151 #endif
1330 jgs 94 binaryOp(right,plus<double>());
1331     return (*this);
1332     }
1333    
1334 jgs 102 Data&
1335     Data::operator-=(const Data& right)
1336 jgs 94 {
1337 jgs 151 #if defined DOPROF
1338 jgs 123 profData->binary++;
1339 jgs 151 #endif
1340 jgs 94 binaryOp(right,minus<double>());
1341     return (*this);
1342     }
1343    
1344 jgs 102 Data&
1345     Data::operator-=(const boost::python::object& right)
1346 jgs 94 {
1347 jgs 151 #if defined DOPROF
1348 jgs 123 profData->binary++;
1349 jgs 151 #endif
1350 jgs 94 binaryOp(right,minus<double>());
1351     return (*this);
1352     }
1353    
1354 jgs 102 Data&
1355     Data::operator*=(const Data& right)
1356 jgs 94 {
1357 jgs 151 #if defined DOPROF
1358 jgs 123 profData->binary++;
1359 jgs 151 #endif
1360 jgs 94 binaryOp(right,multiplies<double>());
1361     return (*this);
1362     }
1363    
1364 jgs 102 Data&
1365     Data::operator*=(const boost::python::object& right)
1366 jgs 94 {
1367 jgs 151 #if defined DOPROF
1368 jgs 123 profData->binary++;
1369 jgs 151 #endif
1370 jgs 94 binaryOp(right,multiplies<double>());
1371     return (*this);
1372     }
1373    
1374 jgs 102 Data&
1375     Data::operator/=(const Data& right)
1376 jgs 94 {
1377 jgs 151 #if defined DOPROF
1378 jgs 123 profData->binary++;
1379 jgs 151 #endif
1380 jgs 94 binaryOp(right,divides<double>());
1381     return (*this);
1382     }
1383    
1384 jgs 102 Data&
1385     Data::operator/=(const boost::python::object& right)
1386 jgs 94 {
1387 jgs 151 #if defined DOPROF
1388 jgs 123 profData->binary++;
1389 jgs 151 #endif
1390 jgs 94 binaryOp(right,divides<double>());
1391     return (*this);
1392     }
1393    
1394 jgs 102 Data
1395 gross 699 Data::rpowO(const boost::python::object& left) const
1396     {
1397     #if defined DOPROF
1398     profData->binary++;
1399     #endif
1400     Data left_d(left,*this);
1401     return left_d.powD(*this);
1402     }
1403    
1404     Data
1405 jgs 102 Data::powO(const boost::python::object& right) const
1406 jgs 94 {
1407 jgs 151 #if defined DOPROF
1408 jgs 123 profData->binary++;
1409 jgs 151 #endif
1410 jgs 94 Data result;
1411     result.copy(*this);
1412     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1413     return result;
1414     }
1415    
1416 jgs 102 Data
1417     Data::powD(const Data& right) const
1418 jgs 94 {
1419 jgs 151 #if defined DOPROF
1420 jgs 123 profData->binary++;
1421 jgs 151 #endif
1422 jgs 94 Data result;
1423     result.copy(*this);
1424     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1425     return result;
1426     }
1427    
1428     //
1429 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1430 jgs 102 Data
1431     escript::operator+(const Data& left, const Data& right)
1432 jgs 94 {
1433     Data result;
1434     //
1435     // perform a deep copy
1436     result.copy(left);
1437     result+=right;
1438     return result;
1439     }
1440    
1441     //
1442 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1443 jgs 102 Data
1444     escript::operator-(const Data& left, const Data& right)
1445 jgs 94 {
1446     Data result;
1447     //
1448     // perform a deep copy
1449     result.copy(left);
1450     result-=right;
1451     return result;
1452     }
1453    
1454     //
1455 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1456 jgs 102 Data
1457     escript::operator*(const Data& left, const Data& right)
1458 jgs 94 {
1459     Data result;
1460     //
1461     // perform a deep copy
1462     result.copy(left);
1463     result*=right;
1464     return result;
1465     }
1466    
1467     //
1468 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1469 jgs 102 Data
1470     escript::operator/(const Data& left, const Data& right)
1471 jgs 94 {
1472     Data result;
1473     //
1474     // perform a deep copy
1475     result.copy(left);
1476     result/=right;
1477     return result;
1478     }
1479    
1480     //
1481 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1482 jgs 102 Data
1483     escript::operator+(const Data& left, const boost::python::object& right)
1484 jgs 94 {
1485     //
1486     // Convert to DataArray format if possible
1487     DataArray temp(right);
1488     Data result;
1489     //
1490     // perform a deep copy
1491     result.copy(left);
1492     result+=right;
1493     return result;
1494     }
1495    
1496     //
1497 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1498 jgs 102 Data
1499     escript::operator-(const Data& left, const boost::python::object& right)
1500 jgs 94 {
1501     //
1502     // Convert to DataArray format if possible
1503     DataArray temp(right);
1504     Data result;
1505     //
1506     // perform a deep copy
1507     result.copy(left);
1508     result-=right;
1509     return result;
1510     }
1511    
1512     //
1513 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1514 jgs 102 Data
1515     escript::operator*(const Data& left, const boost::python::object& right)
1516 jgs 94 {
1517     //
1518     // Convert to DataArray format if possible
1519     DataArray temp(right);
1520     Data result;
1521     //
1522     // perform a deep copy
1523     result.copy(left);
1524     result*=right;
1525     return result;
1526     }
1527    
1528     //
1529 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1530 jgs 102 Data
1531     escript::operator/(const Data& left, const boost::python::object& right)
1532 jgs 94 {
1533     //
1534     // Convert to DataArray format if possible
1535     DataArray temp(right);
1536     Data result;
1537     //
1538     // perform a deep copy
1539     result.copy(left);
1540     result/=right;
1541     return result;
1542     }
1543    
1544     //
1545 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1546 jgs 102 Data
1547     escript::operator+(const boost::python::object& left, const Data& right)
1548 jgs 94 {
1549     //
1550     // Construct the result using the given value and the other parameters
1551     // from right
1552     Data result(left,right);
1553     result+=right;
1554     return result;
1555     }
1556    
1557     //
1558 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1559 jgs 102 Data
1560     escript::operator-(const boost::python::object& left, const Data& right)
1561 jgs 94 {
1562     //
1563     // Construct the result using the given value and the other parameters
1564     // from right
1565     Data result(left,right);
1566     result-=right;
1567     return result;
1568     }
1569    
1570     //
1571 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1572 jgs 102 Data
1573     escript::operator*(const boost::python::object& left, const Data& right)
1574 jgs 94 {
1575     //
1576     // Construct the result using the given value and the other parameters
1577     // from right
1578     Data result(left,right);
1579     result*=right;
1580     return result;
1581     }
1582    
1583     //
1584 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1585 jgs 102 Data
1586     escript::operator/(const boost::python::object& left, const Data& right)
1587 jgs 94 {
1588     //
1589     // Construct the result using the given value and the other parameters
1590     // from right
1591     Data result(left,right);
1592     result/=right;
1593     return result;
1594     }
1595    
1596     //
1597 jgs 102 //bool escript::operator==(const Data& left, const Data& right)
1598     //{
1599     // /*
1600     // NB: this operator does very little at this point, and isn't to
1601     // be relied on. Requires further implementation.
1602     // */
1603     //
1604     // bool ret;
1605     //
1606     // if (left.isEmpty()) {
1607     // if(!right.isEmpty()) {
1608     // ret = false;
1609     // } else {
1610     // ret = true;
1611     // }
1612     // }
1613     //
1614     // if (left.isConstant()) {
1615     // if(!right.isConstant()) {
1616     // ret = false;
1617     // } else {
1618     // ret = true;
1619     // }
1620     // }
1621     //
1622     // if (left.isTagged()) {
1623     // if(!right.isTagged()) {
1624     // ret = false;
1625     // } else {
1626     // ret = true;
1627     // }
1628     // }
1629     //
1630     // if (left.isExpanded()) {
1631     // if(!right.isExpanded()) {
1632     // ret = false;
1633     // } else {
1634     // ret = true;
1635     // }
1636     // }
1637     //
1638     // return ret;
1639     //}
1640    
1641     Data
1642     Data::getItem(const boost::python::object& key) const
1643 jgs 94 {
1644 jgs 102 const DataArrayView& view=getPointDataView();
1645 jgs 94
1646 jgs 102 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1647 jgs 94
1648 jgs 102 if (slice_region.size()!=view.getRank()) {
1649     throw DataException("Error - slice size does not match Data rank.");
1650 jgs 94 }
1651    
1652 jgs 102 return getSlice(slice_region);
1653 jgs 94 }
1654    
1655     Data
1656 jgs 102 Data::getSlice(const DataArrayView::RegionType& region) const
1657 jgs 94 {
1658 jgs 151 #if defined DOPROF
1659 jgs 123 profData->slicing++;
1660 jgs 151 #endif
1661 jgs 102 return Data(*this,region);
1662 jgs 94 }
1663    
1664     void
1665 jgs 102 Data::setItemO(const boost::python::object& key,
1666     const boost::python::object& value)
1667 jgs 94 {
1668 jgs 102 Data tempData(value,getFunctionSpace());
1669     setItemD(key,tempData);
1670     }
1671    
1672     void
1673     Data::setItemD(const boost::python::object& key,
1674     const Data& value)
1675     {
1676 jgs 94 const DataArrayView& view=getPointDataView();
1677 jgs 123
1678 jgs 94 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1679     if (slice_region.size()!=view.getRank()) {
1680     throw DataException("Error - slice size does not match Data rank.");
1681     }
1682 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
1683     setSlice(Data(value,getFunctionSpace()),slice_region);
1684     } else {
1685     setSlice(value,slice_region);
1686     }
1687 jgs 94 }
1688    
1689     void
1690 jgs 102 Data::setSlice(const Data& value,
1691     const DataArrayView::RegionType& region)
1692 jgs 94 {
1693 jgs 151 #if defined DOPROF
1694 jgs 123 profData->slicing++;
1695 jgs 151 #endif
1696 jgs 102 Data tempValue(value);
1697     typeMatchLeft(tempValue);
1698     typeMatchRight(tempValue);
1699     m_data->setSlice(tempValue.m_data.get(),region);
1700     }
1701    
1702     void
1703     Data::typeMatchLeft(Data& right) const
1704     {
1705     if (isExpanded()){
1706     right.expand();
1707     } else if (isTagged()) {
1708     if (right.isConstant()) {
1709     right.tag();
1710     }
1711     }
1712     }
1713    
1714     void
1715     Data::typeMatchRight(const Data& right)
1716     {
1717 jgs 94 if (isTagged()) {
1718     if (right.isExpanded()) {
1719     expand();
1720     }
1721     } else if (isConstant()) {
1722     if (right.isExpanded()) {
1723     expand();
1724     } else if (right.isTagged()) {
1725     tag();
1726     }
1727     }
1728     }
1729    
1730     void
1731     Data::setTaggedValue(int tagKey,
1732     const boost::python::object& value)
1733     {
1734     //
1735     // Ensure underlying data object is of type DataTagged
1736     tag();
1737    
1738     if (!isTagged()) {
1739     throw DataException("Error - DataTagged conversion failed!!");
1740     }
1741    
1742     //
1743     // Construct DataArray from boost::python::object input value
1744     DataArray valueDataArray(value);
1745    
1746     //
1747     // Call DataAbstract::setTaggedValue
1748     m_data->setTaggedValue(tagKey,valueDataArray.getView());
1749     }
1750    
1751 jgs 110 void
1752 jgs 121 Data::setTaggedValueFromCPP(int tagKey,
1753     const DataArrayView& value)
1754     {
1755     //
1756     // Ensure underlying data object is of type DataTagged
1757     tag();
1758    
1759     if (!isTagged()) {
1760     throw DataException("Error - DataTagged conversion failed!!");
1761     }
1762    
1763     //
1764     // Call DataAbstract::setTaggedValue
1765     m_data->setTaggedValue(tagKey,value);
1766     }
1767    
1768 jgs 149 int
1769     Data::getTagNumber(int dpno)
1770     {
1771     return m_data->getTagNumber(dpno);
1772     }
1773    
1774 jgs 121 void
1775 jgs 110 Data::setRefValue(int ref,
1776     const boost::python::numeric::array& value)
1777     {
1778     //
1779     // Construct DataArray from boost::python::object input value
1780     DataArray valueDataArray(value);
1781    
1782     //
1783     // Call DataAbstract::setRefValue
1784     m_data->setRefValue(ref,valueDataArray);
1785     }
1786    
1787     void
1788     Data::getRefValue(int ref,
1789     boost::python::numeric::array& value)
1790     {
1791     //
1792     // Construct DataArray for boost::python::object return value
1793     DataArray valueDataArray(value);
1794    
1795     //
1796     // Load DataArray with values from data-points specified by ref
1797     m_data->getRefValue(ref,valueDataArray);
1798    
1799     //
1800     // Load values from valueDataArray into return numarray
1801    
1802     // extract the shape of the numarray
1803     int rank = value.getrank();
1804     DataArrayView::ShapeType shape;
1805     for (int i=0; i < rank; i++) {
1806     shape.push_back(extract<int>(value.getshape()[i]));
1807     }
1808    
1809     // and load the numarray with the data from the DataArray
1810     DataArrayView valueView = valueDataArray.getView();
1811    
1812     if (rank==0) {
1813 jgs 126 boost::python::numeric::array temp_numArray(valueView());
1814     value = temp_numArray;
1815 jgs 110 }
1816     if (rank==1) {
1817     for (int i=0; i < shape[0]; i++) {
1818     value[i] = valueView(i);
1819     }
1820     }
1821     if (rank==2) {
1822 jgs 126 for (int i=0; i < shape[0]; i++) {
1823     for (int j=0; j < shape[1]; j++) {
1824     value[i][j] = valueView(i,j);
1825     }
1826     }
1827 jgs 110 }
1828     if (rank==3) {
1829 jgs 126 for (int i=0; i < shape[0]; i++) {
1830     for (int j=0; j < shape[1]; j++) {
1831     for (int k=0; k < shape[2]; k++) {
1832     value[i][j][k] = valueView(i,j,k);
1833     }
1834     }
1835     }
1836 jgs 110 }
1837     if (rank==4) {
1838 jgs 126 for (int i=0; i < shape[0]; i++) {
1839     for (int j=0; j < shape[1]; j++) {
1840     for (int k=0; k < shape[2]; k++) {
1841     for (int l=0; l < shape[3]; l++) {
1842     value[i][j][k][l] = valueView(i,j,k,l);
1843     }
1844     }
1845     }
1846     }
1847 jgs 110 }
1848    
1849     }
1850    
1851 jgs 94 void
1852 jgs 119 Data::archiveData(const std::string fileName)
1853     {
1854     cout << "Archiving Data object to: " << fileName << endl;
1855    
1856     //
1857     // Determine type of this Data object
1858     int dataType = -1;
1859    
1860     if (isEmpty()) {
1861     dataType = 0;
1862     cout << "\tdataType: DataEmpty" << endl;
1863     }
1864     if (isConstant()) {
1865     dataType = 1;
1866     cout << "\tdataType: DataConstant" << endl;
1867     }
1868     if (isTagged()) {
1869     dataType = 2;
1870     cout << "\tdataType: DataTagged" << endl;
1871     }
1872     if (isExpanded()) {
1873     dataType = 3;
1874     cout << "\tdataType: DataExpanded" << endl;
1875     }
1876 jgs 123
1877 jgs 119 if (dataType == -1) {
1878     throw DataException("archiveData Error: undefined dataType");
1879     }
1880    
1881     //
1882     // Collect data items common to all Data types
1883     int noSamples = getNumSamples();
1884     int noDPPSample = getNumDataPointsPerSample();
1885     int functionSpaceType = getFunctionSpace().getTypeCode();
1886     int dataPointRank = getDataPointRank();
1887     int dataPointSize = getDataPointSize();
1888     int dataLength = getLength();
1889     DataArrayView::ShapeType dataPointShape = getDataPointShape();
1890     int referenceNumbers[noSamples];
1891     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1892     referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1893     }
1894     int tagNumbers[noSamples];
1895     if (isTagged()) {
1896     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1897     tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1898     }
1899     }
1900    
1901     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1902     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1903     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1904    
1905     //
1906     // Flatten Shape to an array of integers suitable for writing to file
1907     int flatShape[4] = {0,0,0,0};
1908     cout << "\tshape: < ";
1909     for (int dim=0; dim<dataPointRank; dim++) {
1910     flatShape[dim] = dataPointShape[dim];
1911     cout << dataPointShape[dim] << " ";
1912     }
1913     cout << ">" << endl;
1914    
1915     //
1916 jgs 123 // Open archive file
1917 jgs 119 ofstream archiveFile;
1918     archiveFile.open(fileName.data(), ios::out);
1919    
1920     if (!archiveFile.good()) {
1921     throw DataException("archiveData Error: problem opening archive file");
1922     }
1923    
1924 jgs 123 //
1925     // Write common data items to archive file
1926 jgs 119 archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1927     archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1928     archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1929     archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1930     archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1931     archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1932     archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1933     for (int dim = 0; dim < 4; dim++) {
1934     archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1935     }
1936     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1937     archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1938     }
1939     if (isTagged()) {
1940     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1941     archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1942     }
1943     }
1944    
1945     if (!archiveFile.good()) {
1946     throw DataException("archiveData Error: problem writing to archive file");
1947     }
1948    
1949     //
1950 jgs 123 // Archive underlying data values for each Data type
1951     int noValues;
1952 jgs 119 switch (dataType) {
1953     case 0:
1954     // DataEmpty
1955 jgs 123 noValues = 0;
1956     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1957     cout << "\tnoValues: " << noValues << endl;
1958 jgs 119 break;
1959     case 1:
1960     // DataConstant
1961 jgs 123 noValues = m_data->getLength();
1962     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1963     cout << "\tnoValues: " << noValues << endl;
1964     if (m_data->archiveData(archiveFile,noValues)) {
1965     throw DataException("archiveData Error: problem writing data to archive file");
1966     }
1967 jgs 119 break;
1968     case 2:
1969     // DataTagged
1970 jgs 123 noValues = m_data->getLength();
1971     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1972     cout << "\tnoValues: " << noValues << endl;
1973     if (m_data->archiveData(archiveFile,noValues)) {
1974     throw DataException("archiveData Error: problem writing data to archive file");
1975     }
1976 jgs 119 break;
1977     case 3:
1978     // DataExpanded
1979 jgs 123 noValues = m_data->getLength();
1980     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1981     cout << "\tnoValues: " << noValues << endl;
1982     if (m_data->archiveData(archiveFile,noValues)) {
1983     throw DataException("archiveData Error: problem writing data to archive file");
1984     }
1985 jgs 119 break;
1986     }
1987    
1988 jgs 123 if (!archiveFile.good()) {
1989     throw DataException("archiveData Error: problem writing data to archive file");
1990     }
1991    
1992     //
1993     // Close archive file
1994     archiveFile.close();
1995    
1996     if (!archiveFile.good()) {
1997     throw DataException("archiveData Error: problem closing archive file");
1998     }
1999    
2000 jgs 119 }
2001    
2002     void
2003     Data::extractData(const std::string fileName,
2004     const FunctionSpace& fspace)
2005     {
2006     //
2007     // Can only extract Data to an object which is initially DataEmpty
2008     if (!isEmpty()) {
2009     throw DataException("extractData Error: can only extract to DataEmpty object");
2010     }
2011    
2012     cout << "Extracting Data object from: " << fileName << endl;
2013    
2014     int dataType;
2015     int noSamples;
2016     int noDPPSample;
2017     int functionSpaceType;
2018     int dataPointRank;
2019     int dataPointSize;
2020     int dataLength;
2021     DataArrayView::ShapeType dataPointShape;
2022     int flatShape[4];
2023    
2024     //
2025 jgs 123 // Open the archive file
2026 jgs 119 ifstream archiveFile;
2027     archiveFile.open(fileName.data(), ios::in);
2028    
2029     if (!archiveFile.good()) {
2030     throw DataException("extractData Error: problem opening archive file");
2031     }
2032    
2033 jgs 123 //
2034     // Read common data items from archive file
2035 jgs 119 archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
2036     archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2037     archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2038     archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2039     archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2040     archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2041     archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
2042     for (int dim = 0; dim < 4; dim++) {
2043     archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2044     if (flatShape[dim]>0) {
2045     dataPointShape.push_back(flatShape[dim]);
2046     }
2047     }
2048     int referenceNumbers[noSamples];
2049     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2050     archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2051     }
2052     int tagNumbers[noSamples];
2053     if (dataType==2) {
2054     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2055     archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2056     }
2057     }
2058    
2059     if (!archiveFile.good()) {
2060     throw DataException("extractData Error: problem reading from archive file");
2061     }
2062    
2063 jgs 123 //
2064     // Verify the values just read from the archive file
2065 jgs 119 switch (dataType) {
2066     case 0:
2067     cout << "\tdataType: DataEmpty" << endl;
2068     break;
2069     case 1:
2070     cout << "\tdataType: DataConstant" << endl;
2071     break;
2072     case 2:
2073     cout << "\tdataType: DataTagged" << endl;
2074     break;
2075     case 3:
2076     cout << "\tdataType: DataExpanded" << endl;
2077     break;
2078     default:
2079     throw DataException("extractData Error: undefined dataType read from archive file");
2080     break;
2081     }
2082    
2083     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2084     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2085     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2086     cout << "\tshape: < ";
2087     for (int dim = 0; dim < dataPointRank; dim++) {
2088     cout << dataPointShape[dim] << " ";
2089     }
2090     cout << ">" << endl;
2091    
2092     //
2093     // Verify that supplied FunctionSpace object is compatible with this Data object.
2094     if ( (fspace.getTypeCode()!=functionSpaceType) ||
2095     (fspace.getNumSamples()!=noSamples) ||
2096     (fspace.getNumDPPSample()!=noDPPSample)
2097     ) {
2098     throw DataException("extractData Error: incompatible FunctionSpace");
2099     }
2100     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2101     if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2102     throw DataException("extractData Error: incompatible FunctionSpace");
2103     }
2104     }
2105     if (dataType==2) {
2106     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2107     if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2108     throw DataException("extractData Error: incompatible FunctionSpace");
2109     }
2110     }
2111     }
2112    
2113     //
2114     // Construct a DataVector to hold underlying data values
2115     DataVector dataVec(dataLength);
2116    
2117     //
2118     // Load this DataVector with the appropriate values
2119 jgs 123 int noValues;
2120     archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2121     cout << "\tnoValues: " << noValues << endl;
2122 jgs 119 switch (dataType) {
2123     case 0:
2124     // DataEmpty
2125 jgs 123 if (noValues != 0) {
2126     throw DataException("extractData Error: problem reading data from archive file");
2127     }
2128 jgs 119 break;
2129     case 1:
2130     // DataConstant
2131 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2132     throw DataException("extractData Error: problem reading data from archive file");
2133     }
2134 jgs 119 break;
2135     case 2:
2136     // DataTagged
2137 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2138     throw DataException("extractData Error: problem reading data from archive file");
2139     }
2140 jgs 119 break;
2141     case 3:
2142     // DataExpanded
2143 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2144     throw DataException("extractData Error: problem reading data from archive file");
2145     }
2146 jgs 119 break;
2147     }
2148    
2149 jgs 123 if (!archiveFile.good()) {
2150     throw DataException("extractData Error: problem reading from archive file");
2151     }
2152    
2153 jgs 119 //
2154 jgs 123 // Close archive file
2155     archiveFile.close();
2156    
2157     if (!archiveFile.good()) {
2158     throw DataException("extractData Error: problem closing archive file");
2159     }
2160    
2161     //
2162 jgs 119 // Construct an appropriate Data object
2163     DataAbstract* tempData;
2164     switch (dataType) {
2165     case 0:
2166     // DataEmpty
2167     tempData=new DataEmpty();
2168     break;
2169     case 1:
2170     // DataConstant
2171     tempData=new DataConstant(fspace,dataPointShape,dataVec);
2172     break;
2173     case 2:
2174     // DataTagged
2175     tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2176     break;
2177     case 3:
2178     // DataExpanded
2179     tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2180     break;
2181     }
2182     shared_ptr<DataAbstract> temp_data(tempData);
2183     m_data=temp_data;
2184     }
2185    
2186 jgs 94 ostream& escript::operator<<(ostream& o, const Data& data)
2187     {
2188     o << data.toString();
2189     return o;
2190     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26