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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 147 - (hide annotations)
Fri Aug 12 01:45:47 2005 UTC (14 years, 9 months ago) by jgs
File size: 51355 byte(s)
erge of development branch dev-02 back to main trunk on 2005-08-12

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 126 void
481     Data::fillFromNumArray(const boost::python::numeric::array num_array)
482     {
483     //
484     // check rank:
485     if (num_array.getrank()<getDataPointRank())
486     throw DataException("Rank of numarray does not match Data object rank");
487     //
488     // check rank of num_array
489     for (int i=0; i<getDataPointRank(); i++) {
490     if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])
491     throw DataException("Shape of numarray does not match Data object rank");
492     }
493     //
494     // make sure data is expanded:
495     if (! isExpanded()) expand();
496     //
497     // and copy over:
498     m_data->copyAll(num_array);
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     // dataPointView()=numArray[dataPoint];
540     // }
541     // if (dataPointRank==1) {
542     // for (int i=0; i<dataPointShape[0]; i++) {
543     // dataPointView(i)=numArray[dataPoint][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     // }
575     // }
576    
577     //
578     // return the loaded array
579     // return numArray;
580    
581     }
582    
583 jgs 121 const
584 jgs 94 boost::python::numeric::array
585 jgs 117 Data::convertToNumArray()
586     {
587     //
588 jgs 121 // determine the total number of data points
589 jgs 117 int numSamples = getNumSamples();
590     int numDataPointsPerSample = getNumDataPointsPerSample();
591     int numDataPoints = numSamples * numDataPointsPerSample;
592    
593     //
594     // determine the rank and shape of each data point
595     int dataPointRank = getDataPointRank();
596     DataArrayView::ShapeType dataPointShape = getDataPointShape();
597    
598     //
599     // create the numeric array to be returned
600     boost::python::numeric::array numArray(0.0);
601    
602     //
603     // the rank of the returned numeric array will be the rank of
604     // the data points, plus one. Where the rank of the array is n,
605     // the last n-1 dimensions will be equal to the shape of the
606     // data points, whilst the first dimension will be equal to the
607     // total number of data points. Thus the array will consist of
608     // a serial vector of the data points.
609     int arrayRank = dataPointRank + 1;
610     DataArrayView::ShapeType arrayShape;
611     arrayShape.push_back(numDataPoints);
612     for (int d=0; d<dataPointRank; d++) {
613     arrayShape.push_back(dataPointShape[d]);
614     }
615    
616     //
617     // resize the numeric array to the shape just calculated
618     if (arrayRank==1) {
619     numArray.resize(arrayShape[0]);
620     }
621     if (arrayRank==2) {
622     numArray.resize(arrayShape[0],arrayShape[1]);
623     }
624     if (arrayRank==3) {
625     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
626     }
627     if (arrayRank==4) {
628     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
629     }
630     if (arrayRank==5) {
631     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
632     }
633    
634     //
635     // loop through each data point in turn, loading the values for that data point
636     // into the numeric array.
637     int dataPoint = 0;
638     for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
639     for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
640     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
641     if (dataPointRank==0) {
642     numArray[dataPoint]=dataPointView();
643     }
644     if (dataPointRank==1) {
645     for (int i=0; i<dataPointShape[0]; i++) {
646     numArray[dataPoint][i]=dataPointView(i);
647     }
648     }
649     if (dataPointRank==2) {
650     for (int i=0; i<dataPointShape[0]; i++) {
651     for (int j=0; j<dataPointShape[1]; j++) {
652     numArray[dataPoint][i][j] = dataPointView(i,j);
653     }
654     }
655     }
656     if (dataPointRank==3) {
657     for (int i=0; i<dataPointShape[0]; i++) {
658     for (int j=0; j<dataPointShape[1]; j++) {
659     for (int k=0; k<dataPointShape[2]; k++) {
660     numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
661     }
662     }
663     }
664     }
665     if (dataPointRank==4) {
666     for (int i=0; i<dataPointShape[0]; i++) {
667     for (int j=0; j<dataPointShape[1]; j++) {
668     for (int k=0; k<dataPointShape[2]; k++) {
669     for (int l=0; l<dataPointShape[3]; l++) {
670     numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
671     }
672     }
673     }
674     }
675     }
676     dataPoint++;
677 jgs 121 }
678     }
679 jgs 117
680 jgs 121 //
681     // return the loaded array
682     return numArray;
683     }
684    
685     const
686     boost::python::numeric::array
687     Data::convertToNumArrayFromSampleNo(int sampleNo)
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     // determine the number of data points per sample
697     int numDataPointsPerSample = getNumDataPointsPerSample();
698    
699     //
700     // determine the rank and shape of each data point
701     int dataPointRank = getDataPointRank();
702     DataArrayView::ShapeType dataPointShape = getDataPointShape();
703    
704     //
705     // create the numeric array to be returned
706     boost::python::numeric::array numArray(0.0);
707    
708     //
709     // the rank of the returned numeric array will be the rank of
710     // the data points, plus one. Where the rank of the array is n,
711     // the last n-1 dimensions will be equal to the shape of the
712     // data points, whilst the first dimension will be equal to the
713     // total number of data points. Thus the array will consist of
714     // a serial vector of the data points.
715     int arrayRank = dataPointRank + 1;
716     DataArrayView::ShapeType arrayShape;
717     arrayShape.push_back(numDataPointsPerSample);
718     for (int d=0; d<dataPointRank; d++) {
719     arrayShape.push_back(dataPointShape[d]);
720     }
721    
722     //
723     // resize the numeric array to the shape just calculated
724     if (arrayRank==1) {
725     numArray.resize(arrayShape[0]);
726     }
727     if (arrayRank==2) {
728     numArray.resize(arrayShape[0],arrayShape[1]);
729     }
730     if (arrayRank==3) {
731     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
732     }
733     if (arrayRank==4) {
734     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
735     }
736     if (arrayRank==5) {
737     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
738     }
739    
740     //
741     // loop through each data point in turn, loading the values for that data point
742     // into the numeric array.
743     for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {
744     DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);
745     if (dataPointRank==0) {
746     numArray[dataPoint]=dataPointView();
747 jgs 117 }
748 jgs 121 if (dataPointRank==1) {
749     for (int i=0; i<dataPointShape[0]; i++) {
750     numArray[dataPoint][i]=dataPointView(i);
751     }
752     }
753     if (dataPointRank==2) {
754     for (int i=0; i<dataPointShape[0]; i++) {
755     for (int j=0; j<dataPointShape[1]; j++) {
756     numArray[dataPoint][i][j] = dataPointView(i,j);
757     }
758     }
759     }
760     if (dataPointRank==3) {
761     for (int i=0; i<dataPointShape[0]; i++) {
762     for (int j=0; j<dataPointShape[1]; j++) {
763     for (int k=0; k<dataPointShape[2]; k++) {
764     numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
765     }
766     }
767     }
768     }
769     if (dataPointRank==4) {
770     for (int i=0; i<dataPointShape[0]; i++) {
771     for (int j=0; j<dataPointShape[1]; j++) {
772     for (int k=0; k<dataPointShape[2]; k++) {
773     for (int l=0; l<dataPointShape[3]; l++) {
774     numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
775     }
776     }
777     }
778     }
779     }
780 jgs 117 }
781    
782     //
783     // return the loaded array
784     return numArray;
785     }
786    
787 jgs 121 const
788 jgs 117 boost::python::numeric::array
789 jgs 121 Data::convertToNumArrayFromDPNo(int sampleNo,
790     int dataPointNo)
791     {
792     //
793     // Check a valid sample number has been supplied
794     if (sampleNo >= getNumSamples()) {
795     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
796     }
797    
798     //
799     // Check a valid data point number has been supplied
800     if (dataPointNo >= getNumDataPointsPerSample()) {
801     throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");
802     }
803    
804     //
805     // determine the rank and shape of each data point
806     int dataPointRank = getDataPointRank();
807     DataArrayView::ShapeType dataPointShape = getDataPointShape();
808    
809     //
810     // create the numeric array to be returned
811     boost::python::numeric::array numArray(0.0);
812    
813     //
814     // the shape of the returned numeric array will be the same
815     // as that of the data point
816     int arrayRank = dataPointRank;
817     DataArrayView::ShapeType arrayShape = dataPointShape;
818    
819     //
820     // resize the numeric array to the shape just calculated
821     if (arrayRank==0) {
822     numArray.resize(1);
823     }
824     if (arrayRank==1) {
825     numArray.resize(arrayShape[0]);
826     }
827     if (arrayRank==2) {
828     numArray.resize(arrayShape[0],arrayShape[1]);
829     }
830     if (arrayRank==3) {
831     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
832     }
833     if (arrayRank==4) {
834     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
835     }
836    
837     //
838     // load the values for the data point into the numeric array.
839     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
840     if (dataPointRank==0) {
841     numArray[0]=dataPointView();
842     }
843     if (dataPointRank==1) {
844     for (int i=0; i<dataPointShape[0]; i++) {
845     numArray[i]=dataPointView(i);
846     }
847     }
848     if (dataPointRank==2) {
849     for (int i=0; i<dataPointShape[0]; i++) {
850     for (int j=0; j<dataPointShape[1]; j++) {
851     numArray[i][j] = dataPointView(i,j);
852     }
853     }
854     }
855     if (dataPointRank==3) {
856     for (int i=0; i<dataPointShape[0]; i++) {
857     for (int j=0; j<dataPointShape[1]; j++) {
858     for (int k=0; k<dataPointShape[2]; k++) {
859     numArray[i][j][k]=dataPointView(i,j,k);
860     }
861     }
862     }
863     }
864     if (dataPointRank==4) {
865     for (int i=0; i<dataPointShape[0]; i++) {
866     for (int j=0; j<dataPointShape[1]; j++) {
867     for (int k=0; k<dataPointShape[2]; k++) {
868     for (int l=0; l<dataPointShape[3]; l++) {
869     numArray[i][j][k][l]=dataPointView(i,j,k,l);
870     }
871     }
872     }
873     }
874     }
875    
876     //
877     // return the loaded array
878     return numArray;
879     }
880    
881     boost::python::numeric::array
882 jgs 94 Data::integrate() const
883     {
884     int index;
885     int rank = getDataPointRank();
886     DataArrayView::ShapeType shape = getDataPointShape();
887    
888 jgs 123 profData->integrate++;
889    
890 jgs 94 //
891     // calculate the integral values
892     vector<double> integrals(getDataPointSize());
893     AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);
894    
895     //
896     // create the numeric array to be returned
897     // and load the array with the integral values
898     boost::python::numeric::array bp_array(1.0);
899     if (rank==0) {
900 jgs 108 bp_array.resize(1);
901 jgs 94 index = 0;
902     bp_array[0] = integrals[index];
903     }
904     if (rank==1) {
905     bp_array.resize(shape[0]);
906     for (int i=0; i<shape[0]; i++) {
907     index = i;
908     bp_array[i] = integrals[index];
909     }
910     }
911     if (rank==2) {
912     bp_array.resize(shape[0],shape[1]);
913     for (int i=0; i<shape[0]; i++) {
914     for (int j=0; j<shape[1]; j++) {
915     index = i + shape[0] * j;
916     bp_array[i,j] = integrals[index];
917     }
918     }
919     }
920     if (rank==3) {
921     bp_array.resize(shape[0],shape[1],shape[2]);
922     for (int i=0; i<shape[0]; i++) {
923     for (int j=0; j<shape[1]; j++) {
924     for (int k=0; k<shape[2]; k++) {
925     index = i + shape[0] * ( j + shape[1] * k );
926     bp_array[i,j,k] = integrals[index];
927     }
928     }
929     }
930     }
931     if (rank==4) {
932     bp_array.resize(shape[0],shape[1],shape[2],shape[3]);
933     for (int i=0; i<shape[0]; i++) {
934     for (int j=0; j<shape[1]; j++) {
935     for (int k=0; k<shape[2]; k++) {
936     for (int l=0; l<shape[3]; l++) {
937     index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
938     bp_array[i,j,k,l] = integrals[index];
939     }
940     }
941     }
942     }
943     }
944    
945     //
946     // return the loaded array
947     return bp_array;
948     }
949    
950     Data
951     Data::sin() const
952     {
953 jgs 123 profData->unary++;
954 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
955     }
956    
957     Data
958     Data::cos() const
959     {
960 jgs 123 profData->unary++;
961 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
962     }
963    
964     Data
965     Data::tan() const
966     {
967 jgs 123 profData->unary++;
968 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
969     }
970    
971     Data
972     Data::log() const
973     {
974 jgs 123 profData->unary++;
975 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
976     }
977    
978     Data
979     Data::ln() const
980     {
981 jgs 123 profData->unary++;
982 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
983     }
984    
985 jgs 106 Data
986     Data::sign() const
987 jgs 94 {
988 jgs 123 profData->unary++;
989 jgs 106 return escript::unaryOp(*this,escript::fsign);
990 jgs 94 }
991    
992 jgs 106 Data
993     Data::abs() const
994 jgs 94 {
995 jgs 123 profData->unary++;
996 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
997 jgs 94 }
998    
999 jgs 106 Data
1000     Data::neg() const
1001 jgs 94 {
1002 jgs 123 profData->unary++;
1003 jgs 106 return escript::unaryOp(*this,negate<double>());
1004 jgs 94 }
1005    
1006 jgs 102 Data
1007 jgs 106 Data::pos() const
1008 jgs 94 {
1009 jgs 123 profData->unary++;
1010 jgs 102 return (*this);
1011     }
1012    
1013     Data
1014 jgs 106 Data::exp() const
1015 jgs 102 {
1016 jgs 123 profData->unary++;
1017 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1018 jgs 102 }
1019    
1020     Data
1021 jgs 106 Data::sqrt() const
1022 jgs 102 {
1023 jgs 123 profData->unary++;
1024 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1025 jgs 102 }
1026    
1027 jgs 106 double
1028     Data::Lsup() const
1029 jgs 102 {
1030 jgs 123 profData->reduction1++;
1031 jgs 106 //
1032     // set the initial absolute maximum value to zero
1033 jgs 147 AbsMax abs_max_func;
1034     return algorithm(abs_max_func,0);
1035 jgs 102 }
1036    
1037 jgs 106 double
1038 jgs 117 Data::Linf() const
1039     {
1040 jgs 123 profData->reduction1++;
1041 jgs 117 //
1042     // set the initial absolute minimum value to max double
1043 jgs 147 AbsMin abs_min_func;
1044     return algorithm(abs_min_func,numeric_limits<double>::max());
1045 jgs 117 }
1046    
1047     double
1048 jgs 106 Data::sup() const
1049 jgs 102 {
1050 jgs 123 profData->reduction1++;
1051 jgs 106 //
1052     // set the initial maximum value to min possible double
1053 jgs 147 FMax fmax_func;
1054     return algorithm(fmax_func,numeric_limits<double>::max()*-1);
1055 jgs 102 }
1056    
1057 jgs 106 double
1058     Data::inf() const
1059 jgs 102 {
1060 jgs 123 profData->reduction1++;
1061 jgs 106 //
1062     // set the initial minimum value to max possible double
1063 jgs 147 FMin fmin_func;
1064     return algorithm(fmin_func,numeric_limits<double>::max());
1065 jgs 102 }
1066    
1067     Data
1068 jgs 106 Data::maxval() const
1069 jgs 102 {
1070 jgs 123 profData->reduction2++;
1071 jgs 113 //
1072     // set the initial maximum value to min possible double
1073 jgs 147 FMax fmax_func;
1074     return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1075 jgs 102 }
1076    
1077     Data
1078 jgs 106 Data::minval() const
1079 jgs 102 {
1080 jgs 123 profData->reduction2++;
1081 jgs 113 //
1082     // set the initial minimum value to max possible double
1083 jgs 147 FMin fmin_func;
1084     return dp_algorithm(fmin_func,numeric_limits<double>::max());
1085 jgs 102 }
1086    
1087 jgs 123 Data
1088     Data::length() const
1089     {
1090     profData->reduction2++;
1091 jgs 147 Length len_func;
1092     return dp_algorithm(len_func,0);
1093 jgs 123 }
1094    
1095     Data
1096     Data::trace() const
1097     {
1098     profData->reduction2++;
1099 jgs 147 Trace trace_func;
1100     return dp_algorithm(trace_func,0);
1101 jgs 123 }
1102    
1103     Data
1104     Data::transpose(int axis) const
1105     {
1106     profData->reduction2++;
1107     // not implemented
1108     throw DataException("Error - Data::transpose not implemented yet.");
1109     return Data();
1110     }
1111    
1112 jgs 121 const boost::python::tuple
1113     Data::mindp() const
1114     {
1115     Data temp=minval();
1116    
1117     int numSamples=temp.getNumSamples();
1118     int numDPPSample=temp.getNumDataPointsPerSample();
1119    
1120     int i,j,lowi=0,lowj=0;
1121     double min=numeric_limits<double>::max();
1122    
1123     for (i=0; i<numSamples; i++) {
1124     for (j=0; j<numDPPSample; j++) {
1125     double next=temp.getDataPoint(i,j)();
1126     if (next<min) {
1127     min=next;
1128     lowi=i;
1129     lowj=j;
1130     }
1131     }
1132     }
1133    
1134     return make_tuple(lowi,lowj);
1135     }
1136    
1137 jgs 104 void
1138     Data::saveDX(std::string fileName) const
1139     {
1140     getDomain().saveDX(fileName,*this);
1141     return;
1142     }
1143    
1144 jgs 110 void
1145     Data::saveVTK(std::string fileName) const
1146     {
1147     getDomain().saveVTK(fileName,*this);
1148     return;
1149     }
1150    
1151 jgs 102 Data&
1152     Data::operator+=(const Data& right)
1153     {
1154 jgs 123 profData->binary++;
1155 jgs 94 binaryOp(right,plus<double>());
1156     return (*this);
1157     }
1158    
1159 jgs 102 Data&
1160     Data::operator+=(const boost::python::object& right)
1161 jgs 94 {
1162 jgs 123 profData->binary++;
1163 jgs 94 binaryOp(right,plus<double>());
1164     return (*this);
1165     }
1166    
1167 jgs 102 Data&
1168     Data::operator-=(const Data& right)
1169 jgs 94 {
1170 jgs 123 profData->binary++;
1171 jgs 94 binaryOp(right,minus<double>());
1172     return (*this);
1173     }
1174    
1175 jgs 102 Data&
1176     Data::operator-=(const boost::python::object& right)
1177 jgs 94 {
1178 jgs 123 profData->binary++;
1179 jgs 94 binaryOp(right,minus<double>());
1180     return (*this);
1181     }
1182    
1183 jgs 102 Data&
1184     Data::operator*=(const Data& right)
1185 jgs 94 {
1186 jgs 123 profData->binary++;
1187 jgs 94 binaryOp(right,multiplies<double>());
1188     return (*this);
1189     }
1190    
1191 jgs 102 Data&
1192     Data::operator*=(const boost::python::object& right)
1193 jgs 94 {
1194 jgs 123 profData->binary++;
1195 jgs 94 binaryOp(right,multiplies<double>());
1196     return (*this);
1197     }
1198    
1199 jgs 102 Data&
1200     Data::operator/=(const Data& right)
1201 jgs 94 {
1202 jgs 123 profData->binary++;
1203 jgs 94 binaryOp(right,divides<double>());
1204     return (*this);
1205     }
1206    
1207 jgs 102 Data&
1208     Data::operator/=(const boost::python::object& right)
1209 jgs 94 {
1210 jgs 123 profData->binary++;
1211 jgs 94 binaryOp(right,divides<double>());
1212     return (*this);
1213     }
1214    
1215 jgs 102 Data
1216     Data::powO(const boost::python::object& right) const
1217 jgs 94 {
1218 jgs 123 profData->binary++;
1219 jgs 94 Data result;
1220     result.copy(*this);
1221     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1222     return result;
1223     }
1224    
1225 jgs 102 Data
1226     Data::powD(const Data& right) const
1227 jgs 94 {
1228 jgs 123 profData->binary++;
1229 jgs 94 Data result;
1230     result.copy(*this);
1231     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1232     return result;
1233     }
1234    
1235     //
1236 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1237 jgs 102 Data
1238     escript::operator+(const Data& left, const Data& right)
1239 jgs 94 {
1240     Data result;
1241     //
1242     // perform a deep copy
1243     result.copy(left);
1244     result+=right;
1245     return result;
1246     }
1247    
1248     //
1249 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1250 jgs 102 Data
1251     escript::operator-(const Data& left, const Data& right)
1252 jgs 94 {
1253     Data result;
1254     //
1255     // perform a deep copy
1256     result.copy(left);
1257     result-=right;
1258     return result;
1259     }
1260    
1261     //
1262 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1263 jgs 102 Data
1264     escript::operator*(const Data& left, const Data& right)
1265 jgs 94 {
1266     Data result;
1267     //
1268     // perform a deep copy
1269     result.copy(left);
1270     result*=right;
1271     return result;
1272     }
1273    
1274     //
1275 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1276 jgs 102 Data
1277     escript::operator/(const Data& left, const Data& right)
1278 jgs 94 {
1279     Data result;
1280     //
1281     // perform a deep copy
1282     result.copy(left);
1283     result/=right;
1284     return result;
1285     }
1286    
1287     //
1288 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1289 jgs 102 Data
1290     escript::operator+(const Data& left, const boost::python::object& right)
1291 jgs 94 {
1292     //
1293     // Convert to DataArray format if possible
1294     DataArray temp(right);
1295     Data result;
1296     //
1297     // perform a deep copy
1298     result.copy(left);
1299     result+=right;
1300     return result;
1301     }
1302    
1303     //
1304 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1305 jgs 102 Data
1306     escript::operator-(const Data& left, const boost::python::object& right)
1307 jgs 94 {
1308     //
1309     // Convert to DataArray format if possible
1310     DataArray temp(right);
1311     Data result;
1312     //
1313     // perform a deep copy
1314     result.copy(left);
1315     result-=right;
1316     return result;
1317     }
1318    
1319     //
1320 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1321 jgs 102 Data
1322     escript::operator*(const Data& left, const boost::python::object& right)
1323 jgs 94 {
1324     //
1325     // Convert to DataArray format if possible
1326     DataArray temp(right);
1327     Data result;
1328     //
1329     // perform a deep copy
1330     result.copy(left);
1331     result*=right;
1332     return result;
1333     }
1334    
1335     //
1336 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1337 jgs 102 Data
1338     escript::operator/(const Data& left, const boost::python::object& right)
1339 jgs 94 {
1340     //
1341     // Convert to DataArray format if possible
1342     DataArray temp(right);
1343     Data result;
1344     //
1345     // perform a deep copy
1346     result.copy(left);
1347     result/=right;
1348     return result;
1349     }
1350    
1351     //
1352 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1353 jgs 102 Data
1354     escript::operator+(const boost::python::object& left, const Data& right)
1355 jgs 94 {
1356     //
1357     // Construct the result using the given value and the other parameters
1358     // from right
1359     Data result(left,right);
1360     result+=right;
1361     return result;
1362     }
1363    
1364     //
1365 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1366 jgs 102 Data
1367     escript::operator-(const boost::python::object& left, const Data& right)
1368 jgs 94 {
1369     //
1370     // Construct the result using the given value and the other parameters
1371     // from right
1372     Data result(left,right);
1373     result-=right;
1374     return result;
1375     }
1376    
1377     //
1378 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1379 jgs 102 Data
1380     escript::operator*(const boost::python::object& left, const Data& right)
1381 jgs 94 {
1382     //
1383     // Construct the result using the given value and the other parameters
1384     // from right
1385     Data result(left,right);
1386     result*=right;
1387     return result;
1388     }
1389    
1390     //
1391 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1392 jgs 102 Data
1393     escript::operator/(const boost::python::object& left, const Data& right)
1394 jgs 94 {
1395     //
1396     // Construct the result using the given value and the other parameters
1397     // from right
1398     Data result(left,right);
1399     result/=right;
1400     return result;
1401     }
1402    
1403     //
1404 jgs 102 //bool escript::operator==(const Data& left, const Data& right)
1405     //{
1406     // /*
1407     // NB: this operator does very little at this point, and isn't to
1408     // be relied on. Requires further implementation.
1409     // */
1410     //
1411     // bool ret;
1412     //
1413     // if (left.isEmpty()) {
1414     // if(!right.isEmpty()) {
1415     // ret = false;
1416     // } else {
1417     // ret = true;
1418     // }
1419     // }
1420     //
1421     // if (left.isConstant()) {
1422     // if(!right.isConstant()) {
1423     // ret = false;
1424     // } else {
1425     // ret = true;
1426     // }
1427     // }
1428     //
1429     // if (left.isTagged()) {
1430     // if(!right.isTagged()) {
1431     // ret = false;
1432     // } else {
1433     // ret = true;
1434     // }
1435     // }
1436     //
1437     // if (left.isExpanded()) {
1438     // if(!right.isExpanded()) {
1439     // ret = false;
1440     // } else {
1441     // ret = true;
1442     // }
1443     // }
1444     //
1445     // return ret;
1446     //}
1447    
1448     Data
1449     Data::getItem(const boost::python::object& key) const
1450 jgs 94 {
1451 jgs 102 const DataArrayView& view=getPointDataView();
1452 jgs 94
1453 jgs 102 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1454 jgs 94
1455 jgs 102 if (slice_region.size()!=view.getRank()) {
1456     throw DataException("Error - slice size does not match Data rank.");
1457 jgs 94 }
1458    
1459 jgs 102 return getSlice(slice_region);
1460 jgs 94 }
1461    
1462     Data
1463 jgs 102 Data::getSlice(const DataArrayView::RegionType& region) const
1464 jgs 94 {
1465 jgs 123 profData->slicing++;
1466 jgs 102 return Data(*this,region);
1467 jgs 94 }
1468    
1469     void
1470 jgs 102 Data::setItemO(const boost::python::object& key,
1471     const boost::python::object& value)
1472 jgs 94 {
1473 jgs 102 Data tempData(value,getFunctionSpace());
1474     setItemD(key,tempData);
1475     }
1476    
1477     void
1478     Data::setItemD(const boost::python::object& key,
1479     const Data& value)
1480     {
1481 jgs 94 const DataArrayView& view=getPointDataView();
1482 jgs 123
1483 jgs 94 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1484     if (slice_region.size()!=view.getRank()) {
1485     throw DataException("Error - slice size does not match Data rank.");
1486     }
1487 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
1488     setSlice(Data(value,getFunctionSpace()),slice_region);
1489     } else {
1490     setSlice(value,slice_region);
1491     }
1492 jgs 94 }
1493    
1494     void
1495 jgs 102 Data::setSlice(const Data& value,
1496     const DataArrayView::RegionType& region)
1497 jgs 94 {
1498 jgs 123 profData->slicing++;
1499 jgs 102 Data tempValue(value);
1500     typeMatchLeft(tempValue);
1501     typeMatchRight(tempValue);
1502     m_data->setSlice(tempValue.m_data.get(),region);
1503     }
1504    
1505     void
1506     Data::typeMatchLeft(Data& right) const
1507     {
1508     if (isExpanded()){
1509     right.expand();
1510     } else if (isTagged()) {
1511     if (right.isConstant()) {
1512     right.tag();
1513     }
1514     }
1515     }
1516    
1517     void
1518     Data::typeMatchRight(const Data& right)
1519     {
1520 jgs 94 if (isTagged()) {
1521     if (right.isExpanded()) {
1522     expand();
1523     }
1524     } else if (isConstant()) {
1525     if (right.isExpanded()) {
1526     expand();
1527     } else if (right.isTagged()) {
1528     tag();
1529     }
1530     }
1531     }
1532    
1533     void
1534     Data::setTaggedValue(int tagKey,
1535     const boost::python::object& value)
1536     {
1537     //
1538     // Ensure underlying data object is of type DataTagged
1539     tag();
1540    
1541     if (!isTagged()) {
1542     throw DataException("Error - DataTagged conversion failed!!");
1543     }
1544    
1545     //
1546     // Construct DataArray from boost::python::object input value
1547     DataArray valueDataArray(value);
1548    
1549     //
1550     // Call DataAbstract::setTaggedValue
1551     m_data->setTaggedValue(tagKey,valueDataArray.getView());
1552     }
1553    
1554 jgs 110 void
1555 jgs 121 Data::setTaggedValueFromCPP(int tagKey,
1556     const DataArrayView& value)
1557     {
1558     //
1559     // Ensure underlying data object is of type DataTagged
1560     tag();
1561    
1562     if (!isTagged()) {
1563     throw DataException("Error - DataTagged conversion failed!!");
1564     }
1565    
1566     //
1567     // Call DataAbstract::setTaggedValue
1568     m_data->setTaggedValue(tagKey,value);
1569     }
1570    
1571     void
1572 jgs 110 Data::setRefValue(int ref,
1573     const boost::python::numeric::array& value)
1574     {
1575     //
1576     // Construct DataArray from boost::python::object input value
1577     DataArray valueDataArray(value);
1578    
1579     //
1580     // Call DataAbstract::setRefValue
1581     m_data->setRefValue(ref,valueDataArray);
1582     }
1583    
1584     void
1585     Data::getRefValue(int ref,
1586     boost::python::numeric::array& value)
1587     {
1588     //
1589     // Construct DataArray for boost::python::object return value
1590     DataArray valueDataArray(value);
1591    
1592     //
1593     // Load DataArray with values from data-points specified by ref
1594     m_data->getRefValue(ref,valueDataArray);
1595    
1596     //
1597     // Load values from valueDataArray into return numarray
1598    
1599     // extract the shape of the numarray
1600     int rank = value.getrank();
1601     DataArrayView::ShapeType shape;
1602     for (int i=0; i < rank; i++) {
1603     shape.push_back(extract<int>(value.getshape()[i]));
1604     }
1605    
1606     // and load the numarray with the data from the DataArray
1607     DataArrayView valueView = valueDataArray.getView();
1608    
1609     if (rank==0) {
1610 jgs 126 boost::python::numeric::array temp_numArray(valueView());
1611     value = temp_numArray;
1612 jgs 110 }
1613     if (rank==1) {
1614     for (int i=0; i < shape[0]; i++) {
1615     value[i] = valueView(i);
1616     }
1617     }
1618     if (rank==2) {
1619 jgs 126 for (int i=0; i < shape[0]; i++) {
1620     for (int j=0; j < shape[1]; j++) {
1621     value[i][j] = valueView(i,j);
1622     }
1623     }
1624 jgs 110 }
1625     if (rank==3) {
1626 jgs 126 for (int i=0; i < shape[0]; i++) {
1627     for (int j=0; j < shape[1]; j++) {
1628     for (int k=0; k < shape[2]; k++) {
1629     value[i][j][k] = valueView(i,j,k);
1630     }
1631     }
1632     }
1633 jgs 110 }
1634     if (rank==4) {
1635 jgs 126 for (int i=0; i < shape[0]; i++) {
1636     for (int j=0; j < shape[1]; j++) {
1637     for (int k=0; k < shape[2]; k++) {
1638     for (int l=0; l < shape[3]; l++) {
1639     value[i][j][k][l] = valueView(i,j,k,l);
1640     }
1641     }
1642     }
1643     }
1644 jgs 110 }
1645    
1646     }
1647    
1648 jgs 94 void
1649 jgs 119 Data::archiveData(const std::string fileName)
1650     {
1651     cout << "Archiving Data object to: " << fileName << endl;
1652    
1653     //
1654     // Determine type of this Data object
1655     int dataType = -1;
1656    
1657     if (isEmpty()) {
1658     dataType = 0;
1659     cout << "\tdataType: DataEmpty" << endl;
1660     }
1661     if (isConstant()) {
1662     dataType = 1;
1663     cout << "\tdataType: DataConstant" << endl;
1664     }
1665     if (isTagged()) {
1666     dataType = 2;
1667     cout << "\tdataType: DataTagged" << endl;
1668     }
1669     if (isExpanded()) {
1670     dataType = 3;
1671     cout << "\tdataType: DataExpanded" << endl;
1672     }
1673 jgs 123
1674 jgs 119 if (dataType == -1) {
1675     throw DataException("archiveData Error: undefined dataType");
1676     }
1677    
1678     //
1679     // Collect data items common to all Data types
1680     int noSamples = getNumSamples();
1681     int noDPPSample = getNumDataPointsPerSample();
1682     int functionSpaceType = getFunctionSpace().getTypeCode();
1683     int dataPointRank = getDataPointRank();
1684     int dataPointSize = getDataPointSize();
1685     int dataLength = getLength();
1686     DataArrayView::ShapeType dataPointShape = getDataPointShape();
1687     int referenceNumbers[noSamples];
1688     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1689     referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1690     }
1691     int tagNumbers[noSamples];
1692     if (isTagged()) {
1693     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1694     tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1695     }
1696     }
1697    
1698     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1699     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1700     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1701    
1702     //
1703     // Flatten Shape to an array of integers suitable for writing to file
1704     int flatShape[4] = {0,0,0,0};
1705     cout << "\tshape: < ";
1706     for (int dim=0; dim<dataPointRank; dim++) {
1707     flatShape[dim] = dataPointShape[dim];
1708     cout << dataPointShape[dim] << " ";
1709     }
1710     cout << ">" << endl;
1711    
1712     //
1713 jgs 123 // Open archive file
1714 jgs 119 ofstream archiveFile;
1715     archiveFile.open(fileName.data(), ios::out);
1716    
1717     if (!archiveFile.good()) {
1718     throw DataException("archiveData Error: problem opening archive file");
1719     }
1720    
1721 jgs 123 //
1722     // Write common data items to archive file
1723 jgs 119 archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1724     archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1725     archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1726     archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1727     archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1728     archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1729     archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1730     for (int dim = 0; dim < 4; dim++) {
1731     archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1732     }
1733     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1734     archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1735     }
1736     if (isTagged()) {
1737     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1738     archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1739     }
1740     }
1741    
1742     if (!archiveFile.good()) {
1743     throw DataException("archiveData Error: problem writing to archive file");
1744     }
1745    
1746     //
1747 jgs 123 // Archive underlying data values for each Data type
1748     int noValues;
1749 jgs 119 switch (dataType) {
1750     case 0:
1751     // DataEmpty
1752 jgs 123 noValues = 0;
1753     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1754     cout << "\tnoValues: " << noValues << endl;
1755 jgs 119 break;
1756     case 1:
1757     // DataConstant
1758 jgs 123 noValues = m_data->getLength();
1759     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1760     cout << "\tnoValues: " << noValues << endl;
1761     if (m_data->archiveData(archiveFile,noValues)) {
1762     throw DataException("archiveData Error: problem writing data to archive file");
1763     }
1764 jgs 119 break;
1765     case 2:
1766     // DataTagged
1767 jgs 123 noValues = m_data->getLength();
1768     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1769     cout << "\tnoValues: " << noValues << endl;
1770     if (m_data->archiveData(archiveFile,noValues)) {
1771     throw DataException("archiveData Error: problem writing data to archive file");
1772     }
1773 jgs 119 break;
1774     case 3:
1775     // DataExpanded
1776 jgs 123 noValues = m_data->getLength();
1777     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1778     cout << "\tnoValues: " << noValues << endl;
1779     if (m_data->archiveData(archiveFile,noValues)) {
1780     throw DataException("archiveData Error: problem writing data to archive file");
1781     }
1782 jgs 119 break;
1783     }
1784    
1785 jgs 123 if (!archiveFile.good()) {
1786     throw DataException("archiveData Error: problem writing data to archive file");
1787     }
1788    
1789     //
1790     // Close archive file
1791     archiveFile.close();
1792    
1793     if (!archiveFile.good()) {
1794     throw DataException("archiveData Error: problem closing archive file");
1795     }
1796    
1797 jgs 119 }
1798    
1799     void
1800     Data::extractData(const std::string fileName,
1801     const FunctionSpace& fspace)
1802     {
1803     //
1804     // Can only extract Data to an object which is initially DataEmpty
1805     if (!isEmpty()) {
1806     throw DataException("extractData Error: can only extract to DataEmpty object");
1807     }
1808    
1809     cout << "Extracting Data object from: " << fileName << endl;
1810    
1811     int dataType;
1812     int noSamples;
1813     int noDPPSample;
1814     int functionSpaceType;
1815     int dataPointRank;
1816     int dataPointSize;
1817     int dataLength;
1818     DataArrayView::ShapeType dataPointShape;
1819     int flatShape[4];
1820    
1821     //
1822 jgs 123 // Open the archive file
1823 jgs 119 ifstream archiveFile;
1824     archiveFile.open(fileName.data(), ios::in);
1825    
1826     if (!archiveFile.good()) {
1827     throw DataException("extractData Error: problem opening archive file");
1828     }
1829    
1830 jgs 123 //
1831     // Read common data items from archive file
1832 jgs 119 archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1833     archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1834     archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1835     archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1836     archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1837     archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1838     archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1839     for (int dim = 0; dim < 4; dim++) {
1840     archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1841     if (flatShape[dim]>0) {
1842     dataPointShape.push_back(flatShape[dim]);
1843     }
1844     }
1845     int referenceNumbers[noSamples];
1846     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1847     archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1848     }
1849     int tagNumbers[noSamples];
1850     if (dataType==2) {
1851     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1852     archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1853     }
1854     }
1855    
1856     if (!archiveFile.good()) {
1857     throw DataException("extractData Error: problem reading from archive file");
1858     }
1859    
1860 jgs 123 //
1861     // Verify the values just read from the archive file
1862 jgs 119 switch (dataType) {
1863     case 0:
1864     cout << "\tdataType: DataEmpty" << endl;
1865     break;
1866     case 1:
1867     cout << "\tdataType: DataConstant" << endl;
1868     break;
1869     case 2:
1870     cout << "\tdataType: DataTagged" << endl;
1871     break;
1872     case 3:
1873     cout << "\tdataType: DataExpanded" << endl;
1874     break;
1875     default:
1876     throw DataException("extractData Error: undefined dataType read from archive file");
1877     break;
1878     }
1879    
1880     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1881     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1882     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1883     cout << "\tshape: < ";
1884     for (int dim = 0; dim < dataPointRank; dim++) {
1885     cout << dataPointShape[dim] << " ";
1886     }
1887     cout << ">" << endl;
1888    
1889     //
1890     // Verify that supplied FunctionSpace object is compatible with this Data object.
1891     if ( (fspace.getTypeCode()!=functionSpaceType) ||
1892     (fspace.getNumSamples()!=noSamples) ||
1893     (fspace.getNumDPPSample()!=noDPPSample)
1894     ) {
1895     throw DataException("extractData Error: incompatible FunctionSpace");
1896     }
1897     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1898     if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
1899     throw DataException("extractData Error: incompatible FunctionSpace");
1900     }
1901     }
1902     if (dataType==2) {
1903     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1904     if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
1905     throw DataException("extractData Error: incompatible FunctionSpace");
1906     }
1907     }
1908     }
1909    
1910     //
1911     // Construct a DataVector to hold underlying data values
1912     DataVector dataVec(dataLength);
1913    
1914     //
1915     // Load this DataVector with the appropriate values
1916 jgs 123 int noValues;
1917     archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
1918     cout << "\tnoValues: " << noValues << endl;
1919 jgs 119 switch (dataType) {
1920     case 0:
1921     // DataEmpty
1922 jgs 123 if (noValues != 0) {
1923     throw DataException("extractData Error: problem reading data from archive file");
1924     }
1925 jgs 119 break;
1926     case 1:
1927     // DataConstant
1928 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
1929     throw DataException("extractData Error: problem reading data from archive file");
1930     }
1931 jgs 119 break;
1932     case 2:
1933     // DataTagged
1934 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
1935     throw DataException("extractData Error: problem reading data from archive file");
1936     }
1937 jgs 119 break;
1938     case 3:
1939     // DataExpanded
1940 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
1941     throw DataException("extractData Error: problem reading data from archive file");
1942     }
1943 jgs 119 break;
1944     }
1945    
1946 jgs 123 if (!archiveFile.good()) {
1947     throw DataException("extractData Error: problem reading from archive file");
1948     }
1949    
1950 jgs 119 //
1951 jgs 123 // Close archive file
1952     archiveFile.close();
1953    
1954     if (!archiveFile.good()) {
1955     throw DataException("extractData Error: problem closing archive file");
1956     }
1957    
1958     //
1959 jgs 119 // Construct an appropriate Data object
1960     DataAbstract* tempData;
1961     switch (dataType) {
1962     case 0:
1963     // DataEmpty
1964     tempData=new DataEmpty();
1965     break;
1966     case 1:
1967     // DataConstant
1968     tempData=new DataConstant(fspace,dataPointShape,dataVec);
1969     break;
1970     case 2:
1971     // DataTagged
1972     tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
1973     break;
1974     case 3:
1975     // DataExpanded
1976     tempData=new DataExpanded(fspace,dataPointShape,dataVec);
1977     break;
1978     }
1979     shared_ptr<DataAbstract> temp_data(tempData);
1980     m_data=temp_data;
1981     }
1982    
1983 jgs 94 ostream& escript::operator<<(ostream& o, const Data& data)
1984     {
1985     o << data.toString();
1986     return o;
1987     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26