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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26