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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 121 - (hide annotations)
Fri May 6 04:26:16 2005 UTC (14 years, 1 month ago) by jgs
File size: 43372 byte(s)
Merge of development branch back to main trunk on 2005-05-06

1 jgs 94 // $Id$
2     /*=============================================================================
3    
4     ******************************************************************************
5     * *
6     * COPYRIGHT ACcESS 2004 - All Rights Reserved *
7     * *
8     * This software is the property of ACcESS. No part of this code *
9     * may be copied in any form or by any means without the expressed written *
10     * consent of ACcESS. Copying, use or modification of this software *
11     * by any unauthorised person is illegal unless that *
12     * person has a software license agreement with ACcESS. *
13     * *
14     ******************************************************************************
15    
16     ******************************************************************************/
17    
18     #include "escript/Data/Data.h"
19    
20     #include <iostream>
21 jgs 119 #include <fstream>
22 jgs 94 #include <algorithm>
23     #include <vector>
24     #include <exception>
25     #include <functional>
26     #include <math.h>
27    
28     #include <boost/python/str.hpp>
29     #include <boost/python/extract.hpp>
30     #include <boost/python/long.hpp>
31    
32     #include "escript/Data/DataException.h"
33     #include "escript/Data/DataExpanded.h"
34     #include "escript/Data/DataConstant.h"
35     #include "escript/Data/DataTagged.h"
36     #include "escript/Data/DataEmpty.h"
37     #include "escript/Data/DataArray.h"
38     #include "escript/Data/DataAlgorithm.h"
39     #include "escript/Data/FunctionSpaceFactory.h"
40     #include "escript/Data/AbstractContinuousDomain.h"
41 jgs 102 #include "escript/Data/UnaryFuncs.h"
42 jgs 94
43     using namespace std;
44     using namespace boost::python;
45     using namespace boost;
46     using namespace escript;
47    
48     Data::Data()
49     {
50     //
51     // Default data is type DataEmpty
52     DataAbstract* temp=new DataEmpty();
53 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
54     m_data=temp_data;
55 jgs 94 }
56    
57     Data::Data(double value,
58     const tuple& shape,
59     const FunctionSpace& what,
60     bool expanded)
61     {
62     DataArrayView::ShapeType dataPointShape;
63     for (int i = 0; i < shape.attr("__len__")(); ++i) {
64     dataPointShape.push_back(extract<const int>(shape[i]));
65     }
66     DataArray temp(dataPointShape,value);
67     initialise(temp.getView(),what,expanded);
68     }
69    
70     Data::Data(double value,
71     const DataArrayView::ShapeType& dataPointShape,
72     const FunctionSpace& what,
73     bool expanded)
74     {
75     DataArray temp(dataPointShape,value);
76     pair<int,int> dataShape=what.getDataShape();
77     initialise(temp.getView(),what,expanded);
78     }
79    
80 jgs 102 Data::Data(const Data& inData)
81 jgs 94 {
82 jgs 102 m_data=inData.m_data;
83 jgs 94 }
84    
85     Data::Data(const Data& inData,
86     const DataArrayView::RegionType& region)
87     {
88     //
89 jgs 102 // Create Data which is a slice of another Data
90     DataAbstract* tmp = inData.m_data->getSlice(region);
91     shared_ptr<DataAbstract> temp_data(tmp);
92     m_data=temp_data;
93 jgs 94 }
94    
95     Data::Data(const Data& inData,
96     const FunctionSpace& functionspace)
97     {
98     if (inData.getFunctionSpace()==functionspace) {
99     m_data=inData.m_data;
100     } else {
101     Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
102     // Note for Lutz, Must use a reference or pointer to a derived object
103     // in order to get polymorphic behaviour. Shouldn't really
104     // be able to create an instance of AbstractDomain but that was done
105     // as a boost python work around which may no longer be required.
106     const AbstractDomain& inDataDomain=inData.getDomain();
107     if (inDataDomain==functionspace.getDomain()) {
108     inDataDomain.interpolateOnDomain(tmp,inData);
109     } else {
110     inDataDomain.interpolateACross(tmp,inData);
111     }
112     m_data=tmp.m_data;
113     }
114     }
115    
116     Data::Data(const DataTagged::TagListType& tagKeys,
117     const DataTagged::ValueListType & values,
118     const DataArrayView& defaultValue,
119     const FunctionSpace& what,
120     bool expanded)
121     {
122     DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
123 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
124     m_data=temp_data;
125 jgs 94 if (expanded) {
126     expand();
127     }
128     }
129    
130     Data::Data(const numeric::array& value,
131     const FunctionSpace& what,
132     bool expanded)
133     {
134     initialise(value,what,expanded);
135     }
136    
137     Data::Data(const DataArrayView& value,
138     const FunctionSpace& what,
139     bool expanded)
140     {
141     initialise(value,what,expanded);
142     }
143    
144     Data::Data(const object& value,
145     const FunctionSpace& what,
146     bool expanded)
147     {
148     numeric::array asNumArray(value);
149     initialise(asNumArray,what,expanded);
150     }
151    
152     Data::Data(const object& value,
153     const Data& other)
154     {
155     //
156     // Create DataConstant using the given value and all other parameters
157     // copied from other. If value is a rank 0 object this Data
158     // will assume the point data shape of other.
159     DataArray temp(value);
160     if (temp.getView().getRank()==0) {
161     //
162     // Create a DataArray with the scalar value for all elements
163     DataArray temp2(other.getPointDataView().getShape(),temp.getView()());
164     initialise(temp2.getView(),other.getFunctionSpace(),false);
165     } else {
166     //
167     // Create a DataConstant with the same sample shape as other
168     initialise(temp.getView(),other.getFunctionSpace(),false);
169     }
170     }
171    
172     escriptDataC
173     Data::getDataC()
174     {
175     escriptDataC temp;
176     temp.m_dataPtr=(void*)this;
177     return temp;
178     }
179    
180     escriptDataC
181     Data::getDataC() const
182     {
183     escriptDataC temp;
184     temp.m_dataPtr=(void*)this;
185     return temp;
186     }
187    
188 jgs 121 const boost::python::tuple
189 jgs 94 Data::getShapeTuple() const
190     {
191     const DataArrayView::ShapeType& shape=getDataPointShape();
192     switch(getDataPointRank()) {
193     case 0:
194     return make_tuple();
195     case 1:
196     return make_tuple(long_(shape[0]));
197     case 2:
198     return make_tuple(long_(shape[0]),long_(shape[1]));
199     case 3:
200     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]));
201     case 4:
202     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]),long_(shape[3]));
203     default:
204     throw DataException("Error - illegal Data rank.");
205     }
206     }
207    
208     void
209     Data::copy(const Data& other)
210     {
211     //
212     // Perform a deep copy
213     {
214     DataExpanded* temp=dynamic_cast<DataExpanded*>(other.m_data.get());
215     if (temp!=0) {
216     //
217     // Construct a DataExpanded copy
218     DataAbstract* newData=new DataExpanded(*temp);
219 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
220     m_data=temp_data;
221 jgs 94 return;
222     }
223     }
224     {
225     DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());
226     if (temp!=0) {
227     //
228 jgs 102 // Construct a DataTagged copy
229 jgs 94 DataAbstract* newData=new DataTagged(*temp);
230 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
231     m_data=temp_data;
232 jgs 94 return;
233     }
234     }
235     {
236     DataConstant* temp=dynamic_cast<DataConstant*>(other.m_data.get());
237     if (temp!=0) {
238     //
239     // Construct a DataConstant copy
240     DataAbstract* newData=new DataConstant(*temp);
241 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
242     m_data=temp_data;
243 jgs 94 return;
244     }
245     }
246 jgs 102 {
247     DataEmpty* temp=dynamic_cast<DataEmpty*>(other.m_data.get());
248     if (temp!=0) {
249     //
250     // Construct a DataEmpty copy
251     DataAbstract* newData=new DataEmpty();
252     shared_ptr<DataAbstract> temp_data(newData);
253     m_data=temp_data;
254     return;
255     }
256     }
257 jgs 94 throw DataException("Error - Copy not implemented for this Data type.");
258     }
259    
260     void
261     Data::copyWithMask(const Data& other,
262     const Data& mask)
263     {
264     Data mask1;
265     Data mask2;
266    
267     mask1 = mask.wherePositive();
268     mask2.copy(mask1);
269    
270     mask1 *= other;
271     mask2 *= *this;
272     mask2 = *this - mask2;
273    
274     *this = mask1 + mask2;
275     }
276    
277     bool
278     Data::isExpanded() const
279     {
280     DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
281     return (temp!=0);
282     }
283    
284     bool
285     Data::isTagged() const
286     {
287     DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
288     return (temp!=0);
289     }
290    
291     bool
292     Data::isEmpty() const
293     {
294     DataEmpty* temp=dynamic_cast<DataEmpty*>(m_data.get());
295     return (temp!=0);
296     }
297    
298     bool
299     Data::isConstant() const
300     {
301     DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
302     return (temp!=0);
303     }
304    
305     void
306     Data::expand()
307     {
308     if (isConstant()) {
309     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
310     DataAbstract* temp=new DataExpanded(*tempDataConst);
311 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
312     m_data=temp_data;
313 jgs 94 } else if (isTagged()) {
314     DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
315     DataAbstract* temp=new DataExpanded(*tempDataTag);
316 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
317     m_data=temp_data;
318 jgs 94 } else if (isExpanded()) {
319     //
320     // do nothing
321     } else if (isEmpty()) {
322     throw DataException("Error - Expansion of DataEmpty not possible.");
323     } else {
324     throw DataException("Error - Expansion not implemented for this Data type.");
325     }
326     }
327    
328     void
329     Data::tag()
330     {
331     if (isConstant()) {
332     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
333     DataAbstract* temp=new DataTagged(*tempDataConst);
334 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
335     m_data=temp_data;
336 jgs 94 } else if (isTagged()) {
337     // do nothing
338     } else if (isExpanded()) {
339     throw DataException("Error - Creating tag data from DataExpanded not possible.");
340     } else if (isEmpty()) {
341     throw DataException("Error - Creating tag data from DataEmpty not possible.");
342     } else {
343     throw DataException("Error - Tagging not implemented for this Data type.");
344     }
345     }
346    
347 jgs 102 void
348     Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)
349     {
350     m_data->reshapeDataPoint(shape);
351     }
352    
353 jgs 94 Data
354     Data::wherePositive() const
355     {
356     return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
357     }
358    
359     Data
360 jgs 102 Data::whereNegative() const
361     {
362     return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
363     }
364    
365     Data
366 jgs 94 Data::whereNonNegative() const
367     {
368     return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
369     }
370    
371     Data
372 jgs 102 Data::whereNonPositive() const
373 jgs 94 {
374 jgs 102 return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
375 jgs 94 }
376    
377     Data
378     Data::whereZero() const
379     {
380     return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));
381     }
382    
383     Data
384 jgs 102 Data::whereNonZero() const
385     {
386     return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));
387     }
388    
389     Data
390 jgs 94 Data::interpolate(const FunctionSpace& functionspace) const
391     {
392     return Data(*this,functionspace);
393     }
394    
395     bool
396     Data::probeInterpolation(const FunctionSpace& functionspace) const
397     {
398     if (getFunctionSpace()==functionspace) {
399     return true;
400     } else {
401     const AbstractDomain& domain=getDomain();
402     if (domain==functionspace.getDomain()) {
403     return domain.probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());
404     } else {
405     return domain.probeInterpolationACross(getFunctionSpace().getTypeCode(),functionspace.getDomain(),functionspace.getTypeCode());
406     }
407     }
408     }
409    
410     Data
411     Data::gradOn(const FunctionSpace& functionspace) const
412     {
413     if (functionspace.getDomain()!=getDomain())
414     throw DataException("Error - gradient cannot be calculated on different domains.");
415     DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
416     grad_shape.push_back(functionspace.getDim());
417     Data out(0.0,grad_shape,functionspace,true);
418     getDomain().setToGradient(out,*this);
419     return out;
420     }
421    
422     Data
423     Data::grad() const
424     {
425     return gradOn(escript::function(getDomain()));
426     }
427    
428     int
429     Data::getDataPointSize() const
430     {
431     return getPointDataView().noValues();
432     }
433    
434     DataArrayView::ValueType::size_type
435     Data::getLength() const
436     {
437     return m_data->getLength();
438     }
439    
440     const DataArrayView::ShapeType&
441     Data::getDataPointShape() const
442     {
443     return getPointDataView().getShape();
444     }
445    
446 jgs 121 const
447 jgs 94 boost::python::numeric::array
448 jgs 117 Data::convertToNumArray()
449     {
450     //
451 jgs 121 // determine the total number of data points
452 jgs 117 int numSamples = getNumSamples();
453     int numDataPointsPerSample = getNumDataPointsPerSample();
454     int numDataPoints = numSamples * numDataPointsPerSample;
455    
456     //
457     // determine the rank and shape of each data point
458     int dataPointRank = getDataPointRank();
459     DataArrayView::ShapeType dataPointShape = getDataPointShape();
460    
461     //
462     // create the numeric array to be returned
463     boost::python::numeric::array numArray(0.0);
464    
465     //
466     // the rank of the returned numeric array will be the rank of
467     // the data points, plus one. Where the rank of the array is n,
468     // the last n-1 dimensions will be equal to the shape of the
469     // data points, whilst the first dimension will be equal to the
470     // total number of data points. Thus the array will consist of
471     // a serial vector of the data points.
472     int arrayRank = dataPointRank + 1;
473     DataArrayView::ShapeType arrayShape;
474     arrayShape.push_back(numDataPoints);
475     for (int d=0; d<dataPointRank; d++) {
476     arrayShape.push_back(dataPointShape[d]);
477     }
478    
479     //
480     // resize the numeric array to the shape just calculated
481     if (arrayRank==1) {
482     numArray.resize(arrayShape[0]);
483     }
484     if (arrayRank==2) {
485     numArray.resize(arrayShape[0],arrayShape[1]);
486     }
487     if (arrayRank==3) {
488     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
489     }
490     if (arrayRank==4) {
491     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
492     }
493     if (arrayRank==5) {
494     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
495     }
496    
497     //
498     // loop through each data point in turn, loading the values for that data point
499     // into the numeric array.
500     int dataPoint = 0;
501     for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
502     for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
503     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
504     if (dataPointRank==0) {
505     numArray[dataPoint]=dataPointView();
506     }
507     if (dataPointRank==1) {
508     for (int i=0; i<dataPointShape[0]; i++) {
509     numArray[dataPoint][i]=dataPointView(i);
510     }
511     }
512     if (dataPointRank==2) {
513     for (int i=0; i<dataPointShape[0]; i++) {
514     for (int j=0; j<dataPointShape[1]; j++) {
515     numArray[dataPoint][i][j] = dataPointView(i,j);
516     }
517     }
518     }
519     if (dataPointRank==3) {
520     for (int i=0; i<dataPointShape[0]; i++) {
521     for (int j=0; j<dataPointShape[1]; j++) {
522     for (int k=0; k<dataPointShape[2]; k++) {
523     numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
524     }
525     }
526     }
527     }
528     if (dataPointRank==4) {
529     for (int i=0; i<dataPointShape[0]; i++) {
530     for (int j=0; j<dataPointShape[1]; j++) {
531     for (int k=0; k<dataPointShape[2]; k++) {
532     for (int l=0; l<dataPointShape[3]; l++) {
533     numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
534     }
535     }
536     }
537     }
538     }
539     dataPoint++;
540 jgs 121 }
541     }
542 jgs 117
543 jgs 121 //
544     // return the loaded array
545     return numArray;
546     }
547    
548     const
549     boost::python::numeric::array
550     Data::convertToNumArrayFromSampleNo(int sampleNo)
551     {
552     //
553     // Check a valid sample number has been supplied
554     if (sampleNo >= getNumSamples()) {
555     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
556     }
557    
558     //
559     // determine the number of data points per sample
560     int numDataPointsPerSample = getNumDataPointsPerSample();
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(numDataPointsPerSample);
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     for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {
607     DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);
608     if (dataPointRank==0) {
609     numArray[dataPoint]=dataPointView();
610 jgs 117 }
611 jgs 121 if (dataPointRank==1) {
612     for (int i=0; i<dataPointShape[0]; i++) {
613     numArray[dataPoint][i]=dataPointView(i);
614     }
615     }
616     if (dataPointRank==2) {
617     for (int i=0; i<dataPointShape[0]; i++) {
618     for (int j=0; j<dataPointShape[1]; j++) {
619     numArray[dataPoint][i][j] = dataPointView(i,j);
620     }
621     }
622     }
623     if (dataPointRank==3) {
624     for (int i=0; i<dataPointShape[0]; i++) {
625     for (int j=0; j<dataPointShape[1]; j++) {
626     for (int k=0; k<dataPointShape[2]; k++) {
627     numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
628     }
629     }
630     }
631     }
632     if (dataPointRank==4) {
633     for (int i=0; i<dataPointShape[0]; i++) {
634     for (int j=0; j<dataPointShape[1]; j++) {
635     for (int k=0; k<dataPointShape[2]; k++) {
636     for (int l=0; l<dataPointShape[3]; l++) {
637     numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
638     }
639     }
640     }
641     }
642     }
643 jgs 117 }
644    
645     //
646     // return the loaded array
647     return numArray;
648     }
649    
650 jgs 121 const
651 jgs 117 boost::python::numeric::array
652 jgs 121 Data::convertToNumArrayFromDPNo(int sampleNo,
653     int dataPointNo)
654     {
655     //
656     // Check a valid sample number has been supplied
657     if (sampleNo >= getNumSamples()) {
658     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
659     }
660    
661     //
662     // Check a valid data point number has been supplied
663     if (dataPointNo >= getNumDataPointsPerSample()) {
664     throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");
665     }
666    
667     //
668     // determine the rank and shape of each data point
669     int dataPointRank = getDataPointRank();
670     DataArrayView::ShapeType dataPointShape = getDataPointShape();
671    
672     //
673     // create the numeric array to be returned
674     boost::python::numeric::array numArray(0.0);
675    
676     //
677     // the shape of the returned numeric array will be the same
678     // as that of the data point
679     int arrayRank = dataPointRank;
680     DataArrayView::ShapeType arrayShape = dataPointShape;
681    
682     //
683     // resize the numeric array to the shape just calculated
684     if (arrayRank==0) {
685     numArray.resize(1);
686     }
687     if (arrayRank==1) {
688     numArray.resize(arrayShape[0]);
689     }
690     if (arrayRank==2) {
691     numArray.resize(arrayShape[0],arrayShape[1]);
692     }
693     if (arrayRank==3) {
694     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
695     }
696     if (arrayRank==4) {
697     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
698     }
699    
700     //
701     // load the values for the data point into the numeric array.
702     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
703     if (dataPointRank==0) {
704     numArray[0]=dataPointView();
705     }
706     if (dataPointRank==1) {
707     for (int i=0; i<dataPointShape[0]; i++) {
708     numArray[i]=dataPointView(i);
709     }
710     }
711     if (dataPointRank==2) {
712     for (int i=0; i<dataPointShape[0]; i++) {
713     for (int j=0; j<dataPointShape[1]; j++) {
714     numArray[i][j] = dataPointView(i,j);
715     }
716     }
717     }
718     if (dataPointRank==3) {
719     for (int i=0; i<dataPointShape[0]; i++) {
720     for (int j=0; j<dataPointShape[1]; j++) {
721     for (int k=0; k<dataPointShape[2]; k++) {
722     numArray[i][j][k]=dataPointView(i,j,k);
723     }
724     }
725     }
726     }
727     if (dataPointRank==4) {
728     for (int i=0; i<dataPointShape[0]; i++) {
729     for (int j=0; j<dataPointShape[1]; j++) {
730     for (int k=0; k<dataPointShape[2]; k++) {
731     for (int l=0; l<dataPointShape[3]; l++) {
732     numArray[i][j][k][l]=dataPointView(i,j,k,l);
733     }
734     }
735     }
736     }
737     }
738    
739     //
740     // return the loaded array
741     return numArray;
742     }
743    
744     boost::python::numeric::array
745 jgs 94 Data::integrate() const
746     {
747     int index;
748     int rank = getDataPointRank();
749     DataArrayView::ShapeType shape = getDataPointShape();
750    
751     //
752     // calculate the integral values
753     vector<double> integrals(getDataPointSize());
754     AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);
755    
756     //
757     // create the numeric array to be returned
758     // and load the array with the integral values
759     boost::python::numeric::array bp_array(1.0);
760     if (rank==0) {
761 jgs 108 bp_array.resize(1);
762 jgs 94 index = 0;
763     bp_array[0] = integrals[index];
764     }
765     if (rank==1) {
766     bp_array.resize(shape[0]);
767     for (int i=0; i<shape[0]; i++) {
768     index = i;
769     bp_array[i] = integrals[index];
770     }
771     }
772     if (rank==2) {
773     bp_array.resize(shape[0],shape[1]);
774     for (int i=0; i<shape[0]; i++) {
775     for (int j=0; j<shape[1]; j++) {
776     index = i + shape[0] * j;
777     bp_array[i,j] = integrals[index];
778     }
779     }
780     }
781     if (rank==3) {
782     bp_array.resize(shape[0],shape[1],shape[2]);
783     for (int i=0; i<shape[0]; i++) {
784     for (int j=0; j<shape[1]; j++) {
785     for (int k=0; k<shape[2]; k++) {
786     index = i + shape[0] * ( j + shape[1] * k );
787     bp_array[i,j,k] = integrals[index];
788     }
789     }
790     }
791     }
792     if (rank==4) {
793     bp_array.resize(shape[0],shape[1],shape[2],shape[3]);
794     for (int i=0; i<shape[0]; i++) {
795     for (int j=0; j<shape[1]; j++) {
796     for (int k=0; k<shape[2]; k++) {
797     for (int l=0; l<shape[3]; l++) {
798     index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
799     bp_array[i,j,k,l] = integrals[index];
800     }
801     }
802     }
803     }
804     }
805    
806     //
807     // return the loaded array
808     return bp_array;
809     }
810    
811     Data
812     Data::sin() const
813     {
814     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
815     }
816    
817     Data
818     Data::cos() const
819     {
820     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
821     }
822    
823     Data
824     Data::tan() const
825     {
826     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
827     }
828    
829     Data
830     Data::log() const
831     {
832     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
833     }
834    
835     Data
836     Data::ln() const
837     {
838     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
839     }
840    
841 jgs 106 Data
842     Data::sign() const
843 jgs 94 {
844 jgs 106 return escript::unaryOp(*this,escript::fsign);
845 jgs 94 }
846    
847 jgs 106 Data
848     Data::abs() const
849 jgs 94 {
850 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
851 jgs 94 }
852    
853 jgs 106 Data
854     Data::neg() const
855 jgs 94 {
856 jgs 106 return escript::unaryOp(*this,negate<double>());
857 jgs 94 }
858    
859 jgs 102 Data
860 jgs 106 Data::pos() const
861 jgs 94 {
862 jgs 102 return (*this);
863     }
864    
865     Data
866 jgs 106 Data::exp() const
867 jgs 102 {
868 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
869 jgs 102 }
870    
871     Data
872 jgs 106 Data::sqrt() const
873 jgs 102 {
874 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
875 jgs 102 }
876    
877 jgs 106 double
878     Data::Lsup() const
879 jgs 102 {
880 jgs 106 //
881     // set the initial absolute maximum value to zero
882     return algorithm(DataAlgorithmAdapter<AbsMax>(0));
883 jgs 102 }
884    
885 jgs 106 double
886 jgs 117 Data::Linf() const
887     {
888     //
889     // set the initial absolute minimum value to max double
890     return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));
891     }
892    
893     double
894 jgs 106 Data::sup() const
895 jgs 102 {
896 jgs 106 //
897     // set the initial maximum value to min possible double
898 jgs 108 return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
899 jgs 102 }
900    
901 jgs 106 double
902     Data::inf() const
903 jgs 102 {
904 jgs 106 //
905     // set the initial minimum value to max possible double
906     return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));
907 jgs 102 }
908    
909     Data
910 jgs 106 Data::maxval() const
911 jgs 102 {
912 jgs 113 //
913     // set the initial maximum value to min possible double
914 jgs 108 return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
915 jgs 102 }
916    
917     Data
918 jgs 106 Data::minval() const
919 jgs 102 {
920 jgs 113 //
921     // set the initial minimum value to max possible double
922 jgs 106 return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));
923 jgs 102 }
924    
925 jgs 121 const boost::python::tuple
926     Data::mindp() const
927     {
928     Data temp=minval();
929    
930     int numSamples=temp.getNumSamples();
931     int numDPPSample=temp.getNumDataPointsPerSample();
932    
933     int i,j,lowi=0,lowj=0;
934     double min=numeric_limits<double>::max();
935    
936     for (i=0; i<numSamples; i++) {
937     for (j=0; j<numDPPSample; j++) {
938     double next=temp.getDataPoint(i,j)();
939     if (next<min) {
940     min=next;
941     lowi=i;
942     lowj=j;
943     }
944     }
945     }
946    
947     return make_tuple(lowi,lowj);
948     }
949    
950 jgs 102 Data
951 jgs 106 Data::length() const
952 jgs 102 {
953 jgs 106 return dp_algorithm(DataAlgorithmAdapter<Length>(0));
954 jgs 102 }
955    
956     Data
957 jgs 106 Data::trace() const
958 jgs 102 {
959 jgs 106 return dp_algorithm(DataAlgorithmAdapter<Trace>(0));
960 jgs 102 }
961    
962     Data
963 jgs 106 Data::transpose(int axis) const
964 jgs 102 {
965 jgs 106 // not implemented
966     throw DataException("Error - Data::transpose not implemented yet.");
967     return Data();
968 jgs 102 }
969    
970 jgs 104 void
971     Data::saveDX(std::string fileName) const
972     {
973     getDomain().saveDX(fileName,*this);
974     return;
975     }
976    
977 jgs 110 void
978     Data::saveVTK(std::string fileName) const
979     {
980     getDomain().saveVTK(fileName,*this);
981     return;
982     }
983    
984 jgs 102 Data&
985     Data::operator+=(const Data& right)
986     {
987 jgs 94 binaryOp(right,plus<double>());
988     return (*this);
989     }
990    
991 jgs 102 Data&
992     Data::operator+=(const boost::python::object& right)
993 jgs 94 {
994     binaryOp(right,plus<double>());
995     return (*this);
996     }
997    
998 jgs 102 Data&
999     Data::operator-=(const Data& right)
1000 jgs 94 {
1001     binaryOp(right,minus<double>());
1002     return (*this);
1003     }
1004    
1005 jgs 102 Data&
1006     Data::operator-=(const boost::python::object& right)
1007 jgs 94 {
1008     binaryOp(right,minus<double>());
1009     return (*this);
1010     }
1011    
1012 jgs 102 Data&
1013     Data::operator*=(const Data& right)
1014 jgs 94 {
1015     binaryOp(right,multiplies<double>());
1016     return (*this);
1017     }
1018    
1019 jgs 102 Data&
1020     Data::operator*=(const boost::python::object& right)
1021 jgs 94 {
1022     binaryOp(right,multiplies<double>());
1023     return (*this);
1024     }
1025    
1026 jgs 102 Data&
1027     Data::operator/=(const Data& right)
1028 jgs 94 {
1029     binaryOp(right,divides<double>());
1030     return (*this);
1031     }
1032    
1033 jgs 102 Data&
1034     Data::operator/=(const boost::python::object& right)
1035 jgs 94 {
1036     binaryOp(right,divides<double>());
1037     return (*this);
1038     }
1039    
1040 jgs 102 Data
1041     Data::powO(const boost::python::object& right) const
1042 jgs 94 {
1043     Data result;
1044     result.copy(*this);
1045     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1046     return result;
1047     }
1048    
1049 jgs 102 Data
1050     Data::powD(const Data& right) const
1051 jgs 94 {
1052     Data result;
1053     result.copy(*this);
1054     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1055     return result;
1056     }
1057    
1058     //
1059     // NOTE: It is essential to specify the namepsace this operator belongs to
1060 jgs 102 Data
1061     escript::operator+(const Data& left, const Data& right)
1062 jgs 94 {
1063     Data result;
1064     //
1065     // perform a deep copy
1066     result.copy(left);
1067     result+=right;
1068     return result;
1069     }
1070    
1071     //
1072     // NOTE: It is essential to specify the namepsace this operator belongs to
1073 jgs 102 Data
1074     escript::operator-(const Data& left, const Data& right)
1075 jgs 94 {
1076     Data result;
1077     //
1078     // perform a deep copy
1079     result.copy(left);
1080     result-=right;
1081     return result;
1082     }
1083    
1084     //
1085     // NOTE: It is essential to specify the namepsace this operator belongs to
1086 jgs 102 Data
1087     escript::operator*(const Data& left, const Data& right)
1088 jgs 94 {
1089     Data result;
1090     //
1091     // perform a deep copy
1092     result.copy(left);
1093     result*=right;
1094     return result;
1095     }
1096    
1097     //
1098     // NOTE: It is essential to specify the namepsace this operator belongs to
1099 jgs 102 Data
1100     escript::operator/(const Data& left, const Data& right)
1101 jgs 94 {
1102     Data result;
1103     //
1104     // perform a deep copy
1105     result.copy(left);
1106     result/=right;
1107     return result;
1108     }
1109    
1110     //
1111     // NOTE: It is essential to specify the namepsace this operator belongs to
1112 jgs 102 Data
1113     escript::operator+(const Data& left, const boost::python::object& right)
1114 jgs 94 {
1115     //
1116     // Convert to DataArray format if possible
1117     DataArray temp(right);
1118     Data result;
1119     //
1120     // perform a deep copy
1121     result.copy(left);
1122     result+=right;
1123     return result;
1124     }
1125    
1126     //
1127     // NOTE: It is essential to specify the namepsace this operator belongs to
1128 jgs 102 Data
1129     escript::operator-(const Data& left, const boost::python::object& right)
1130 jgs 94 {
1131     //
1132     // Convert to DataArray format if possible
1133     DataArray temp(right);
1134     Data result;
1135     //
1136     // perform a deep copy
1137     result.copy(left);
1138     result-=right;
1139     return result;
1140     }
1141    
1142     //
1143     // NOTE: It is essential to specify the namepsace this operator belongs to
1144 jgs 102 Data
1145     escript::operator*(const Data& left, const boost::python::object& right)
1146 jgs 94 {
1147     //
1148     // Convert to DataArray format if possible
1149     DataArray temp(right);
1150     Data result;
1151     //
1152     // perform a deep copy
1153     result.copy(left);
1154     result*=right;
1155     return result;
1156     }
1157    
1158     //
1159     // NOTE: It is essential to specify the namepsace this operator belongs to
1160 jgs 102 Data
1161     escript::operator/(const Data& left, const boost::python::object& right)
1162 jgs 94 {
1163     //
1164     // Convert to DataArray format if possible
1165     DataArray temp(right);
1166     Data result;
1167     //
1168     // perform a deep copy
1169     result.copy(left);
1170     result/=right;
1171     return result;
1172     }
1173    
1174     //
1175     // NOTE: It is essential to specify the namepsace this operator belongs to
1176 jgs 102 Data
1177     escript::operator+(const boost::python::object& left, const Data& right)
1178 jgs 94 {
1179     //
1180     // Construct the result using the given value and the other parameters
1181     // from right
1182     Data result(left,right);
1183     result+=right;
1184     return result;
1185     }
1186    
1187     //
1188     // NOTE: It is essential to specify the namepsace this operator belongs to
1189 jgs 102 Data
1190     escript::operator-(const boost::python::object& left, const Data& right)
1191 jgs 94 {
1192     //
1193     // Construct the result using the given value and the other parameters
1194     // from right
1195     Data result(left,right);
1196     result-=right;
1197     return result;
1198     }
1199    
1200     //
1201     // NOTE: It is essential to specify the namepsace this operator belongs to
1202 jgs 102 Data
1203     escript::operator*(const boost::python::object& left, const Data& right)
1204 jgs 94 {
1205     //
1206     // Construct the result using the given value and the other parameters
1207     // from right
1208     Data result(left,right);
1209     result*=right;
1210     return result;
1211     }
1212    
1213     //
1214     // NOTE: It is essential to specify the namepsace this operator belongs to
1215 jgs 102 Data
1216     escript::operator/(const boost::python::object& left, const Data& right)
1217 jgs 94 {
1218     //
1219     // Construct the result using the given value and the other parameters
1220     // from right
1221     Data result(left,right);
1222     result/=right;
1223     return result;
1224     }
1225    
1226     //
1227     // NOTE: It is essential to specify the namepsace this operator belongs to
1228 jgs 102 //bool escript::operator==(const Data& left, const Data& right)
1229     //{
1230     // /*
1231     // NB: this operator does very little at this point, and isn't to
1232     // be relied on. Requires further implementation.
1233     // */
1234     //
1235     // bool ret;
1236     //
1237     // if (left.isEmpty()) {
1238     // if(!right.isEmpty()) {
1239     // ret = false;
1240     // } else {
1241     // ret = true;
1242     // }
1243     // }
1244     //
1245     // if (left.isConstant()) {
1246     // if(!right.isConstant()) {
1247     // ret = false;
1248     // } else {
1249     // ret = true;
1250     // }
1251     // }
1252     //
1253     // if (left.isTagged()) {
1254     // if(!right.isTagged()) {
1255     // ret = false;
1256     // } else {
1257     // ret = true;
1258     // }
1259     // }
1260     //
1261     // if (left.isExpanded()) {
1262     // if(!right.isExpanded()) {
1263     // ret = false;
1264     // } else {
1265     // ret = true;
1266     // }
1267     // }
1268     //
1269     // return ret;
1270     //}
1271    
1272     Data
1273     Data::getItem(const boost::python::object& key) const
1274 jgs 94 {
1275 jgs 102 const DataArrayView& view=getPointDataView();
1276 jgs 94
1277 jgs 102 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1278 jgs 94
1279 jgs 102 if (slice_region.size()!=view.getRank()) {
1280     throw DataException("Error - slice size does not match Data rank.");
1281 jgs 94 }
1282    
1283 jgs 102 return getSlice(slice_region);
1284 jgs 94 }
1285    
1286     Data
1287 jgs 102 Data::getSlice(const DataArrayView::RegionType& region) const
1288 jgs 94 {
1289 jgs 102 return Data(*this,region);
1290 jgs 94 }
1291    
1292     void
1293 jgs 102 Data::setItemO(const boost::python::object& key,
1294     const boost::python::object& value)
1295 jgs 94 {
1296 jgs 102 Data tempData(value,getFunctionSpace());
1297     setItemD(key,tempData);
1298     }
1299    
1300     void
1301     Data::setItemD(const boost::python::object& key,
1302     const Data& value)
1303     {
1304 jgs 94 const DataArrayView& view=getPointDataView();
1305     DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1306     if (slice_region.size()!=view.getRank()) {
1307     throw DataException("Error - slice size does not match Data rank.");
1308     }
1309 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
1310     setSlice(Data(value,getFunctionSpace()),slice_region);
1311     } else {
1312     setSlice(value,slice_region);
1313     }
1314 jgs 94 }
1315    
1316     void
1317 jgs 102 Data::setSlice(const Data& value,
1318     const DataArrayView::RegionType& region)
1319 jgs 94 {
1320 jgs 102 Data tempValue(value);
1321     typeMatchLeft(tempValue);
1322     typeMatchRight(tempValue);
1323     m_data->setSlice(tempValue.m_data.get(),region);
1324     }
1325    
1326     void
1327     Data::typeMatchLeft(Data& right) const
1328     {
1329     if (isExpanded()){
1330     right.expand();
1331     } else if (isTagged()) {
1332     if (right.isConstant()) {
1333     right.tag();
1334     }
1335     }
1336     }
1337    
1338     void
1339     Data::typeMatchRight(const Data& right)
1340     {
1341 jgs 94 if (isTagged()) {
1342     if (right.isExpanded()) {
1343     expand();
1344     }
1345     } else if (isConstant()) {
1346     if (right.isExpanded()) {
1347     expand();
1348     } else if (right.isTagged()) {
1349     tag();
1350     }
1351     }
1352     }
1353    
1354     void
1355     Data::setTaggedValue(int tagKey,
1356     const boost::python::object& value)
1357     {
1358     //
1359     // Ensure underlying data object is of type DataTagged
1360     tag();
1361    
1362     if (!isTagged()) {
1363     throw DataException("Error - DataTagged conversion failed!!");
1364     }
1365    
1366     //
1367     // Construct DataArray from boost::python::object input value
1368     DataArray valueDataArray(value);
1369    
1370     //
1371     // Call DataAbstract::setTaggedValue
1372     m_data->setTaggedValue(tagKey,valueDataArray.getView());
1373     }
1374    
1375 jgs 110 void
1376 jgs 121 Data::setTaggedValueFromCPP(int tagKey,
1377     const DataArrayView& value)
1378     {
1379     //
1380     // Ensure underlying data object is of type DataTagged
1381     tag();
1382    
1383     if (!isTagged()) {
1384     throw DataException("Error - DataTagged conversion failed!!");
1385     }
1386    
1387     //
1388     // Call DataAbstract::setTaggedValue
1389     m_data->setTaggedValue(tagKey,value);
1390     }
1391    
1392     void
1393 jgs 110 Data::setRefValue(int ref,
1394     const boost::python::numeric::array& value)
1395     {
1396     //
1397     // Construct DataArray from boost::python::object input value
1398     DataArray valueDataArray(value);
1399    
1400     //
1401     // Call DataAbstract::setRefValue
1402     m_data->setRefValue(ref,valueDataArray);
1403     }
1404    
1405     void
1406     Data::getRefValue(int ref,
1407     boost::python::numeric::array& value)
1408     {
1409     //
1410     // Construct DataArray for boost::python::object return value
1411     DataArray valueDataArray(value);
1412    
1413     //
1414     // Load DataArray with values from data-points specified by ref
1415     m_data->getRefValue(ref,valueDataArray);
1416    
1417     //
1418     // Load values from valueDataArray into return numarray
1419    
1420     // extract the shape of the numarray
1421     int rank = value.getrank();
1422     DataArrayView::ShapeType shape;
1423     for (int i=0; i < rank; i++) {
1424     shape.push_back(extract<int>(value.getshape()[i]));
1425     }
1426    
1427     // and load the numarray with the data from the DataArray
1428     DataArrayView valueView = valueDataArray.getView();
1429    
1430     if (rank==0) {
1431     throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1432     }
1433     if (rank==1) {
1434     for (int i=0; i < shape[0]; i++) {
1435     value[i] = valueView(i);
1436     }
1437     }
1438     if (rank==2) {
1439     throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1440     }
1441     if (rank==3) {
1442     throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1443     }
1444     if (rank==4) {
1445     throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1446     }
1447    
1448     }
1449    
1450 jgs 94 void
1451 jgs 119 Data::archiveData(const std::string fileName)
1452     {
1453     cout << "Archiving Data object to: " << fileName << endl;
1454    
1455     //
1456     // Determine type of this Data object
1457     int dataType = -1;
1458    
1459     if (isEmpty()) {
1460     dataType = 0;
1461     cout << "\tdataType: DataEmpty" << endl;
1462     }
1463     if (isConstant()) {
1464     dataType = 1;
1465     cout << "\tdataType: DataConstant" << endl;
1466     }
1467     if (isTagged()) {
1468     dataType = 2;
1469     cout << "\tdataType: DataTagged" << endl;
1470     }
1471     if (isExpanded()) {
1472     dataType = 3;
1473     cout << "\tdataType: DataExpanded" << endl;
1474     }
1475     if (dataType == -1) {
1476     throw DataException("archiveData Error: undefined dataType");
1477     }
1478    
1479     //
1480     // Collect data items common to all Data types
1481     int noSamples = getNumSamples();
1482     int noDPPSample = getNumDataPointsPerSample();
1483     int functionSpaceType = getFunctionSpace().getTypeCode();
1484     int dataPointRank = getDataPointRank();
1485     int dataPointSize = getDataPointSize();
1486     int dataLength = getLength();
1487     DataArrayView::ShapeType dataPointShape = getDataPointShape();
1488     int referenceNumbers[noSamples];
1489     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1490     referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1491     }
1492     int tagNumbers[noSamples];
1493     if (isTagged()) {
1494     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1495     tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1496     }
1497     }
1498    
1499     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1500     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1501     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1502    
1503     //
1504     // Flatten Shape to an array of integers suitable for writing to file
1505     int flatShape[4] = {0,0,0,0};
1506     cout << "\tshape: < ";
1507     for (int dim=0; dim<dataPointRank; dim++) {
1508     flatShape[dim] = dataPointShape[dim];
1509     cout << dataPointShape[dim] << " ";
1510     }
1511     cout << ">" << endl;
1512    
1513     //
1514     // Write common data items to archive file
1515     ofstream archiveFile;
1516     archiveFile.open(fileName.data(), ios::out);
1517    
1518     if (!archiveFile.good()) {
1519     throw DataException("archiveData Error: problem opening archive file");
1520     }
1521    
1522     archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1523     archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1524     archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1525     archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1526     archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1527     archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1528     archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1529     for (int dim = 0; dim < 4; dim++) {
1530     archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1531     }
1532     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1533     archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1534     }
1535     if (isTagged()) {
1536     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1537     archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1538     }
1539     }
1540    
1541     if (!archiveFile.good()) {
1542     throw DataException("archiveData Error: problem writing to archive file");
1543     }
1544    
1545     archiveFile.close();
1546    
1547     if (!archiveFile.good()) {
1548     throw DataException("archiveData Error: problem closing archive file");
1549     }
1550    
1551     //
1552     // Collect and archive underlying data values for each Data type
1553     switch (dataType) {
1554     case 0:
1555     // DataEmpty
1556     break;
1557     case 1:
1558     // DataConstant
1559     break;
1560     case 2:
1561     // DataTagged
1562     break;
1563     case 3:
1564     // DataExpanded
1565     break;
1566     }
1567    
1568     }
1569    
1570     void
1571     Data::extractData(const std::string fileName,
1572     const FunctionSpace& fspace)
1573     {
1574     //
1575     // Can only extract Data to an object which is initially DataEmpty
1576     if (!isEmpty()) {
1577     throw DataException("extractData Error: can only extract to DataEmpty object");
1578     }
1579    
1580     cout << "Extracting Data object from: " << fileName << endl;
1581    
1582     int dataType;
1583     int noSamples;
1584     int noDPPSample;
1585     int functionSpaceType;
1586     int dataPointRank;
1587     int dataPointSize;
1588     int dataLength;
1589     DataArrayView::ShapeType dataPointShape;
1590     int flatShape[4];
1591    
1592     //
1593     // Open the archive file and read common data items
1594     ifstream archiveFile;
1595     archiveFile.open(fileName.data(), ios::in);
1596    
1597     if (!archiveFile.good()) {
1598     throw DataException("extractData Error: problem opening archive file");
1599     }
1600    
1601     archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1602     archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1603     archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1604     archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1605     archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1606     archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1607     archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1608     for (int dim = 0; dim < 4; dim++) {
1609     archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1610     if (flatShape[dim]>0) {
1611     dataPointShape.push_back(flatShape[dim]);
1612     }
1613     }
1614     int referenceNumbers[noSamples];
1615     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1616     archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1617     }
1618     int tagNumbers[noSamples];
1619     if (dataType==2) {
1620     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1621     archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1622     }
1623     }
1624    
1625     if (!archiveFile.good()) {
1626     throw DataException("extractData Error: problem reading from archive file");
1627     }
1628    
1629     archiveFile.close();
1630    
1631     if (!archiveFile.good()) {
1632     throw DataException("extractData Error: problem closing archive file");
1633     }
1634    
1635     switch (dataType) {
1636     case 0:
1637     cout << "\tdataType: DataEmpty" << endl;
1638     break;
1639     case 1:
1640     cout << "\tdataType: DataConstant" << endl;
1641     break;
1642     case 2:
1643     cout << "\tdataType: DataTagged" << endl;
1644     break;
1645     case 3:
1646     cout << "\tdataType: DataExpanded" << endl;
1647     break;
1648     default:
1649     throw DataException("extractData Error: undefined dataType read from archive file");
1650     break;
1651     }
1652    
1653     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1654     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1655     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1656     cout << "\tshape: < ";
1657     for (int dim = 0; dim < dataPointRank; dim++) {
1658     cout << dataPointShape[dim] << " ";
1659     }
1660     cout << ">" << endl;
1661    
1662     //
1663     // Verify that supplied FunctionSpace object is compatible with this Data object.
1664     if ( (fspace.getTypeCode()!=functionSpaceType) ||
1665     (fspace.getNumSamples()!=noSamples) ||
1666     (fspace.getNumDPPSample()!=noDPPSample)
1667     ) {
1668     throw DataException("extractData Error: incompatible FunctionSpace");
1669     }
1670     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1671     if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
1672     throw DataException("extractData Error: incompatible FunctionSpace");
1673     }
1674     }
1675     if (dataType==2) {
1676     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1677     if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
1678     throw DataException("extractData Error: incompatible FunctionSpace");
1679     }
1680     }
1681     }
1682    
1683     //
1684     // Construct a DataVector to hold underlying data values
1685     DataVector dataVec(dataLength);
1686    
1687     //
1688     // Load this DataVector with the appropriate values
1689     switch (dataType) {
1690     case 0:
1691     // DataEmpty
1692     break;
1693     case 1:
1694     // DataConstant
1695     break;
1696     case 2:
1697     // DataTagged
1698     break;
1699     case 3:
1700     // DataExpanded
1701     break;
1702     }
1703    
1704     //
1705     // Construct an appropriate Data object
1706     DataAbstract* tempData;
1707     switch (dataType) {
1708     case 0:
1709     // DataEmpty
1710     tempData=new DataEmpty();
1711     break;
1712     case 1:
1713     // DataConstant
1714     tempData=new DataConstant(fspace,dataPointShape,dataVec);
1715     break;
1716     case 2:
1717     // DataTagged
1718     tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
1719     break;
1720     case 3:
1721     // DataExpanded
1722     tempData=new DataExpanded(fspace,dataPointShape,dataVec);
1723     break;
1724     }
1725     shared_ptr<DataAbstract> temp_data(tempData);
1726     m_data=temp_data;
1727    
1728     }
1729    
1730 jgs 94 ostream& escript::operator<<(ostream& o, const Data& data)
1731     {
1732     o << data.toString();
1733     return o;
1734     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26