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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 122 - (hide annotations)
Thu Jun 9 05:38:05 2005 UTC (14 years, 3 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/Data.cpp
File size: 43332 byte(s)
Merge of development branch back to main trunk on 2005-06-09

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26