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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26