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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 123 - (hide annotations)
Fri Jul 8 04:08:13 2005 UTC (14 years, 7 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/Data.cpp
File size: 47581 byte(s)
Merge of development branch back to main trunk on 2005-07-08

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26