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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 148 - (hide annotations)
Tue Aug 23 01:24:31 2005 UTC (14 years, 3 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/Data.cpp
File size: 52180 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-08-23

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 148 Data result;
1011     // perform a deep copy
1012     result.copy(*this);
1013     return result;
1014 jgs 102 }
1015    
1016     Data
1017 jgs 106 Data::exp() const
1018 jgs 102 {
1019 jgs 123 profData->unary++;
1020 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1021 jgs 102 }
1022    
1023     Data
1024 jgs 106 Data::sqrt() const
1025 jgs 102 {
1026 jgs 123 profData->unary++;
1027 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1028 jgs 102 }
1029    
1030 jgs 106 double
1031     Data::Lsup() const
1032 jgs 102 {
1033 jgs 123 profData->reduction1++;
1034 jgs 106 //
1035     // set the initial absolute maximum value to zero
1036 jgs 147 AbsMax abs_max_func;
1037     return algorithm(abs_max_func,0);
1038 jgs 102 }
1039    
1040 jgs 106 double
1041 jgs 117 Data::Linf() const
1042     {
1043 jgs 123 profData->reduction1++;
1044 jgs 117 //
1045     // set the initial absolute minimum value to max double
1046 jgs 147 AbsMin abs_min_func;
1047     return algorithm(abs_min_func,numeric_limits<double>::max());
1048 jgs 117 }
1049    
1050     double
1051 jgs 106 Data::sup() const
1052 jgs 102 {
1053 jgs 123 profData->reduction1++;
1054 jgs 106 //
1055     // set the initial maximum value to min possible double
1056 jgs 147 FMax fmax_func;
1057     return algorithm(fmax_func,numeric_limits<double>::max()*-1);
1058 jgs 102 }
1059    
1060 jgs 106 double
1061     Data::inf() const
1062 jgs 102 {
1063 jgs 123 profData->reduction1++;
1064 jgs 106 //
1065     // set the initial minimum value to max possible double
1066 jgs 147 FMin fmin_func;
1067     return algorithm(fmin_func,numeric_limits<double>::max());
1068 jgs 102 }
1069    
1070     Data
1071 jgs 106 Data::maxval() const
1072 jgs 102 {
1073 jgs 123 profData->reduction2++;
1074 jgs 113 //
1075     // set the initial maximum value to min possible double
1076 jgs 147 FMax fmax_func;
1077     return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1078 jgs 102 }
1079    
1080     Data
1081 jgs 106 Data::minval() const
1082 jgs 102 {
1083 jgs 123 profData->reduction2++;
1084 jgs 113 //
1085     // set the initial minimum value to max possible double
1086 jgs 147 FMin fmin_func;
1087     return dp_algorithm(fmin_func,numeric_limits<double>::max());
1088 jgs 102 }
1089    
1090 jgs 123 Data
1091     Data::length() const
1092     {
1093     profData->reduction2++;
1094 jgs 147 Length len_func;
1095     return dp_algorithm(len_func,0);
1096 jgs 123 }
1097    
1098     Data
1099     Data::trace() const
1100     {
1101     profData->reduction2++;
1102 jgs 147 Trace trace_func;
1103     return dp_algorithm(trace_func,0);
1104 jgs 123 }
1105    
1106     Data
1107     Data::transpose(int axis) const
1108     {
1109     profData->reduction2++;
1110     // not implemented
1111     throw DataException("Error - Data::transpose not implemented yet.");
1112     return Data();
1113     }
1114    
1115 jgs 121 const boost::python::tuple
1116     Data::mindp() const
1117     {
1118 jgs 148 // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1119     // abort (for unknown reasons) if there are openmp directives with it in the
1120     // surrounding function
1121    
1122     int SampleNo;
1123     int DataPointNo;
1124    
1125     calc_mindp(SampleNo,DataPointNo);
1126    
1127     return make_tuple(SampleNo,DataPointNo);
1128     }
1129    
1130     void
1131     Data::calc_mindp(int& SampleNo,
1132     int& DataPointNo) const
1133     {
1134     int i,j;
1135     int lowi=0,lowj=0;
1136     double min=numeric_limits<double>::max();
1137    
1138 jgs 121 Data temp=minval();
1139    
1140     int numSamples=temp.getNumSamples();
1141     int numDPPSample=temp.getNumDataPointsPerSample();
1142    
1143 jgs 148 double next,local_min;
1144     int local_lowi,local_lowj;
1145 jgs 121
1146 jgs 148 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1147     {
1148     local_min=min;
1149     #pragma omp for private(i,j) schedule(static)
1150     for (i=0; i<numSamples; i++) {
1151     for (j=0; j<numDPPSample; j++) {
1152     next=temp.getDataPoint(i,j)();
1153     if (next<local_min) {
1154     local_min=next;
1155     local_lowi=i;
1156     local_lowj=j;
1157     }
1158 jgs 121 }
1159     }
1160 jgs 148 #pragma omp critical
1161     if (local_min<min) {
1162     min=local_min;
1163     lowi=local_lowi;
1164     lowj=local_lowj;
1165     }
1166 jgs 121 }
1167    
1168 jgs 148 SampleNo = lowi;
1169     DataPointNo = lowj;
1170 jgs 121 }
1171    
1172 jgs 104 void
1173     Data::saveDX(std::string fileName) const
1174     {
1175     getDomain().saveDX(fileName,*this);
1176     return;
1177     }
1178    
1179 jgs 110 void
1180     Data::saveVTK(std::string fileName) const
1181     {
1182     getDomain().saveVTK(fileName,*this);
1183     return;
1184     }
1185    
1186 jgs 102 Data&
1187     Data::operator+=(const Data& right)
1188     {
1189 jgs 123 profData->binary++;
1190 jgs 94 binaryOp(right,plus<double>());
1191     return (*this);
1192     }
1193    
1194 jgs 102 Data&
1195     Data::operator+=(const boost::python::object& right)
1196 jgs 94 {
1197 jgs 123 profData->binary++;
1198 jgs 94 binaryOp(right,plus<double>());
1199     return (*this);
1200     }
1201    
1202 jgs 102 Data&
1203     Data::operator-=(const Data& right)
1204 jgs 94 {
1205 jgs 123 profData->binary++;
1206 jgs 94 binaryOp(right,minus<double>());
1207     return (*this);
1208     }
1209    
1210 jgs 102 Data&
1211     Data::operator-=(const boost::python::object& right)
1212 jgs 94 {
1213 jgs 123 profData->binary++;
1214 jgs 94 binaryOp(right,minus<double>());
1215     return (*this);
1216     }
1217    
1218 jgs 102 Data&
1219     Data::operator*=(const Data& right)
1220 jgs 94 {
1221 jgs 123 profData->binary++;
1222 jgs 94 binaryOp(right,multiplies<double>());
1223     return (*this);
1224     }
1225    
1226 jgs 102 Data&
1227     Data::operator*=(const boost::python::object& right)
1228 jgs 94 {
1229 jgs 123 profData->binary++;
1230 jgs 94 binaryOp(right,multiplies<double>());
1231     return (*this);
1232     }
1233    
1234 jgs 102 Data&
1235     Data::operator/=(const Data& right)
1236 jgs 94 {
1237 jgs 123 profData->binary++;
1238 jgs 94 binaryOp(right,divides<double>());
1239     return (*this);
1240     }
1241    
1242 jgs 102 Data&
1243     Data::operator/=(const boost::python::object& right)
1244 jgs 94 {
1245 jgs 123 profData->binary++;
1246 jgs 94 binaryOp(right,divides<double>());
1247     return (*this);
1248     }
1249    
1250 jgs 102 Data
1251     Data::powO(const boost::python::object& right) const
1252 jgs 94 {
1253 jgs 123 profData->binary++;
1254 jgs 94 Data result;
1255     result.copy(*this);
1256     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1257     return result;
1258     }
1259    
1260 jgs 102 Data
1261     Data::powD(const Data& right) const
1262 jgs 94 {
1263 jgs 123 profData->binary++;
1264 jgs 94 Data result;
1265     result.copy(*this);
1266     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1267     return result;
1268     }
1269    
1270     //
1271 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1272 jgs 102 Data
1273     escript::operator+(const Data& left, const Data& right)
1274 jgs 94 {
1275     Data result;
1276     //
1277     // perform a deep copy
1278     result.copy(left);
1279     result+=right;
1280     return result;
1281     }
1282    
1283     //
1284 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1285 jgs 102 Data
1286     escript::operator-(const Data& left, const Data& right)
1287 jgs 94 {
1288     Data result;
1289     //
1290     // perform a deep copy
1291     result.copy(left);
1292     result-=right;
1293     return result;
1294     }
1295    
1296     //
1297 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1298 jgs 102 Data
1299     escript::operator*(const Data& left, const Data& right)
1300 jgs 94 {
1301     Data result;
1302     //
1303     // perform a deep copy
1304     result.copy(left);
1305     result*=right;
1306     return result;
1307     }
1308    
1309     //
1310 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1311 jgs 102 Data
1312     escript::operator/(const Data& left, const Data& right)
1313 jgs 94 {
1314     Data result;
1315     //
1316     // perform a deep copy
1317     result.copy(left);
1318     result/=right;
1319     return result;
1320     }
1321    
1322     //
1323 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1324 jgs 102 Data
1325     escript::operator+(const Data& left, const boost::python::object& right)
1326 jgs 94 {
1327     //
1328     // Convert to DataArray format if possible
1329     DataArray temp(right);
1330     Data result;
1331     //
1332     // perform a deep copy
1333     result.copy(left);
1334     result+=right;
1335     return result;
1336     }
1337    
1338     //
1339 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1340 jgs 102 Data
1341     escript::operator-(const Data& left, const boost::python::object& right)
1342 jgs 94 {
1343     //
1344     // Convert to DataArray format if possible
1345     DataArray temp(right);
1346     Data result;
1347     //
1348     // perform a deep copy
1349     result.copy(left);
1350     result-=right;
1351     return result;
1352     }
1353    
1354     //
1355 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1356 jgs 102 Data
1357     escript::operator*(const Data& left, const boost::python::object& right)
1358 jgs 94 {
1359     //
1360     // Convert to DataArray format if possible
1361     DataArray temp(right);
1362     Data result;
1363     //
1364     // perform a deep copy
1365     result.copy(left);
1366     result*=right;
1367     return result;
1368     }
1369    
1370     //
1371 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1372 jgs 102 Data
1373     escript::operator/(const Data& left, const boost::python::object& right)
1374 jgs 94 {
1375     //
1376     // Convert to DataArray format if possible
1377     DataArray temp(right);
1378     Data result;
1379     //
1380     // perform a deep copy
1381     result.copy(left);
1382     result/=right;
1383     return result;
1384     }
1385    
1386     //
1387 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1388 jgs 102 Data
1389     escript::operator+(const boost::python::object& left, const Data& right)
1390 jgs 94 {
1391     //
1392     // Construct the result using the given value and the other parameters
1393     // from right
1394     Data result(left,right);
1395     result+=right;
1396     return result;
1397     }
1398    
1399     //
1400 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1401 jgs 102 Data
1402     escript::operator-(const boost::python::object& left, const Data& right)
1403 jgs 94 {
1404     //
1405     // Construct the result using the given value and the other parameters
1406     // from right
1407     Data result(left,right);
1408     result-=right;
1409     return result;
1410     }
1411    
1412     //
1413 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1414 jgs 102 Data
1415     escript::operator*(const boost::python::object& left, const Data& right)
1416 jgs 94 {
1417     //
1418     // Construct the result using the given value and the other parameters
1419     // from right
1420     Data result(left,right);
1421     result*=right;
1422     return result;
1423     }
1424    
1425     //
1426 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1427 jgs 102 Data
1428     escript::operator/(const boost::python::object& left, const Data& right)
1429 jgs 94 {
1430     //
1431     // Construct the result using the given value and the other parameters
1432     // from right
1433     Data result(left,right);
1434     result/=right;
1435     return result;
1436     }
1437    
1438     //
1439 jgs 102 //bool escript::operator==(const Data& left, const Data& right)
1440     //{
1441     // /*
1442     // NB: this operator does very little at this point, and isn't to
1443     // be relied on. Requires further implementation.
1444     // */
1445     //
1446     // bool ret;
1447     //
1448     // if (left.isEmpty()) {
1449     // if(!right.isEmpty()) {
1450     // ret = false;
1451     // } else {
1452     // ret = true;
1453     // }
1454     // }
1455     //
1456     // if (left.isConstant()) {
1457     // if(!right.isConstant()) {
1458     // ret = false;
1459     // } else {
1460     // ret = true;
1461     // }
1462     // }
1463     //
1464     // if (left.isTagged()) {
1465     // if(!right.isTagged()) {
1466     // ret = false;
1467     // } else {
1468     // ret = true;
1469     // }
1470     // }
1471     //
1472     // if (left.isExpanded()) {
1473     // if(!right.isExpanded()) {
1474     // ret = false;
1475     // } else {
1476     // ret = true;
1477     // }
1478     // }
1479     //
1480     // return ret;
1481     //}
1482    
1483     Data
1484     Data::getItem(const boost::python::object& key) const
1485 jgs 94 {
1486 jgs 102 const DataArrayView& view=getPointDataView();
1487 jgs 94
1488 jgs 102 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1489 jgs 94
1490 jgs 102 if (slice_region.size()!=view.getRank()) {
1491     throw DataException("Error - slice size does not match Data rank.");
1492 jgs 94 }
1493    
1494 jgs 102 return getSlice(slice_region);
1495 jgs 94 }
1496    
1497     Data
1498 jgs 102 Data::getSlice(const DataArrayView::RegionType& region) const
1499 jgs 94 {
1500 jgs 123 profData->slicing++;
1501 jgs 102 return Data(*this,region);
1502 jgs 94 }
1503    
1504     void
1505 jgs 102 Data::setItemO(const boost::python::object& key,
1506     const boost::python::object& value)
1507 jgs 94 {
1508 jgs 102 Data tempData(value,getFunctionSpace());
1509     setItemD(key,tempData);
1510     }
1511    
1512     void
1513     Data::setItemD(const boost::python::object& key,
1514     const Data& value)
1515     {
1516 jgs 94 const DataArrayView& view=getPointDataView();
1517 jgs 123
1518 jgs 94 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1519     if (slice_region.size()!=view.getRank()) {
1520     throw DataException("Error - slice size does not match Data rank.");
1521     }
1522 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
1523     setSlice(Data(value,getFunctionSpace()),slice_region);
1524     } else {
1525     setSlice(value,slice_region);
1526     }
1527 jgs 94 }
1528    
1529     void
1530 jgs 102 Data::setSlice(const Data& value,
1531     const DataArrayView::RegionType& region)
1532 jgs 94 {
1533 jgs 123 profData->slicing++;
1534 jgs 102 Data tempValue(value);
1535     typeMatchLeft(tempValue);
1536     typeMatchRight(tempValue);
1537     m_data->setSlice(tempValue.m_data.get(),region);
1538     }
1539    
1540     void
1541     Data::typeMatchLeft(Data& right) const
1542     {
1543     if (isExpanded()){
1544     right.expand();
1545     } else if (isTagged()) {
1546     if (right.isConstant()) {
1547     right.tag();
1548     }
1549     }
1550     }
1551    
1552     void
1553     Data::typeMatchRight(const Data& right)
1554     {
1555 jgs 94 if (isTagged()) {
1556     if (right.isExpanded()) {
1557     expand();
1558     }
1559     } else if (isConstant()) {
1560     if (right.isExpanded()) {
1561     expand();
1562     } else if (right.isTagged()) {
1563     tag();
1564     }
1565     }
1566     }
1567    
1568     void
1569     Data::setTaggedValue(int tagKey,
1570     const boost::python::object& value)
1571     {
1572     //
1573     // Ensure underlying data object is of type DataTagged
1574     tag();
1575    
1576     if (!isTagged()) {
1577     throw DataException("Error - DataTagged conversion failed!!");
1578     }
1579    
1580     //
1581     // Construct DataArray from boost::python::object input value
1582     DataArray valueDataArray(value);
1583    
1584     //
1585     // Call DataAbstract::setTaggedValue
1586     m_data->setTaggedValue(tagKey,valueDataArray.getView());
1587     }
1588    
1589 jgs 110 void
1590 jgs 121 Data::setTaggedValueFromCPP(int tagKey,
1591     const DataArrayView& value)
1592     {
1593     //
1594     // Ensure underlying data object is of type DataTagged
1595     tag();
1596    
1597     if (!isTagged()) {
1598     throw DataException("Error - DataTagged conversion failed!!");
1599     }
1600    
1601     //
1602     // Call DataAbstract::setTaggedValue
1603     m_data->setTaggedValue(tagKey,value);
1604     }
1605    
1606     void
1607 jgs 110 Data::setRefValue(int ref,
1608     const boost::python::numeric::array& value)
1609     {
1610     //
1611     // Construct DataArray from boost::python::object input value
1612     DataArray valueDataArray(value);
1613    
1614     //
1615     // Call DataAbstract::setRefValue
1616     m_data->setRefValue(ref,valueDataArray);
1617     }
1618    
1619     void
1620     Data::getRefValue(int ref,
1621     boost::python::numeric::array& value)
1622     {
1623     //
1624     // Construct DataArray for boost::python::object return value
1625     DataArray valueDataArray(value);
1626    
1627     //
1628     // Load DataArray with values from data-points specified by ref
1629     m_data->getRefValue(ref,valueDataArray);
1630    
1631     //
1632     // Load values from valueDataArray into return numarray
1633    
1634     // extract the shape of the numarray
1635     int rank = value.getrank();
1636     DataArrayView::ShapeType shape;
1637     for (int i=0; i < rank; i++) {
1638     shape.push_back(extract<int>(value.getshape()[i]));
1639     }
1640    
1641     // and load the numarray with the data from the DataArray
1642     DataArrayView valueView = valueDataArray.getView();
1643    
1644     if (rank==0) {
1645 jgs 126 boost::python::numeric::array temp_numArray(valueView());
1646     value = temp_numArray;
1647 jgs 110 }
1648     if (rank==1) {
1649     for (int i=0; i < shape[0]; i++) {
1650     value[i] = valueView(i);
1651     }
1652     }
1653     if (rank==2) {
1654 jgs 126 for (int i=0; i < shape[0]; i++) {
1655     for (int j=0; j < shape[1]; j++) {
1656     value[i][j] = valueView(i,j);
1657     }
1658     }
1659 jgs 110 }
1660     if (rank==3) {
1661 jgs 126 for (int i=0; i < shape[0]; i++) {
1662     for (int j=0; j < shape[1]; j++) {
1663     for (int k=0; k < shape[2]; k++) {
1664     value[i][j][k] = valueView(i,j,k);
1665     }
1666     }
1667     }
1668 jgs 110 }
1669     if (rank==4) {
1670 jgs 126 for (int i=0; i < shape[0]; i++) {
1671     for (int j=0; j < shape[1]; j++) {
1672     for (int k=0; k < shape[2]; k++) {
1673     for (int l=0; l < shape[3]; l++) {
1674     value[i][j][k][l] = valueView(i,j,k,l);
1675     }
1676     }
1677     }
1678     }
1679 jgs 110 }
1680    
1681     }
1682    
1683 jgs 94 void
1684 jgs 119 Data::archiveData(const std::string fileName)
1685     {
1686     cout << "Archiving Data object to: " << fileName << endl;
1687    
1688     //
1689     // Determine type of this Data object
1690     int dataType = -1;
1691    
1692     if (isEmpty()) {
1693     dataType = 0;
1694     cout << "\tdataType: DataEmpty" << endl;
1695     }
1696     if (isConstant()) {
1697     dataType = 1;
1698     cout << "\tdataType: DataConstant" << endl;
1699     }
1700     if (isTagged()) {
1701     dataType = 2;
1702     cout << "\tdataType: DataTagged" << endl;
1703     }
1704     if (isExpanded()) {
1705     dataType = 3;
1706     cout << "\tdataType: DataExpanded" << endl;
1707     }
1708 jgs 123
1709 jgs 119 if (dataType == -1) {
1710     throw DataException("archiveData Error: undefined dataType");
1711     }
1712    
1713     //
1714     // Collect data items common to all Data types
1715     int noSamples = getNumSamples();
1716     int noDPPSample = getNumDataPointsPerSample();
1717     int functionSpaceType = getFunctionSpace().getTypeCode();
1718     int dataPointRank = getDataPointRank();
1719     int dataPointSize = getDataPointSize();
1720     int dataLength = getLength();
1721     DataArrayView::ShapeType dataPointShape = getDataPointShape();
1722     int referenceNumbers[noSamples];
1723     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1724     referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1725     }
1726     int tagNumbers[noSamples];
1727     if (isTagged()) {
1728     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1729     tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1730     }
1731     }
1732    
1733     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1734     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1735     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1736    
1737     //
1738     // Flatten Shape to an array of integers suitable for writing to file
1739     int flatShape[4] = {0,0,0,0};
1740     cout << "\tshape: < ";
1741     for (int dim=0; dim<dataPointRank; dim++) {
1742     flatShape[dim] = dataPointShape[dim];
1743     cout << dataPointShape[dim] << " ";
1744     }
1745     cout << ">" << endl;
1746    
1747     //
1748 jgs 123 // Open archive file
1749 jgs 119 ofstream archiveFile;
1750     archiveFile.open(fileName.data(), ios::out);
1751    
1752     if (!archiveFile.good()) {
1753     throw DataException("archiveData Error: problem opening archive file");
1754     }
1755    
1756 jgs 123 //
1757     // Write common data items to archive file
1758 jgs 119 archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1759     archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1760     archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1761     archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1762     archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1763     archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1764     archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1765     for (int dim = 0; dim < 4; dim++) {
1766     archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1767     }
1768     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1769     archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1770     }
1771     if (isTagged()) {
1772     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1773     archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1774     }
1775     }
1776    
1777     if (!archiveFile.good()) {
1778     throw DataException("archiveData Error: problem writing to archive file");
1779     }
1780    
1781     //
1782 jgs 123 // Archive underlying data values for each Data type
1783     int noValues;
1784 jgs 119 switch (dataType) {
1785     case 0:
1786     // DataEmpty
1787 jgs 123 noValues = 0;
1788     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1789     cout << "\tnoValues: " << noValues << endl;
1790 jgs 119 break;
1791     case 1:
1792     // DataConstant
1793 jgs 123 noValues = m_data->getLength();
1794     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1795     cout << "\tnoValues: " << noValues << endl;
1796     if (m_data->archiveData(archiveFile,noValues)) {
1797     throw DataException("archiveData Error: problem writing data to archive file");
1798     }
1799 jgs 119 break;
1800     case 2:
1801     // DataTagged
1802 jgs 123 noValues = m_data->getLength();
1803     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1804     cout << "\tnoValues: " << noValues << endl;
1805     if (m_data->archiveData(archiveFile,noValues)) {
1806     throw DataException("archiveData Error: problem writing data to archive file");
1807     }
1808 jgs 119 break;
1809     case 3:
1810     // DataExpanded
1811 jgs 123 noValues = m_data->getLength();
1812     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1813     cout << "\tnoValues: " << noValues << endl;
1814     if (m_data->archiveData(archiveFile,noValues)) {
1815     throw DataException("archiveData Error: problem writing data to archive file");
1816     }
1817 jgs 119 break;
1818     }
1819    
1820 jgs 123 if (!archiveFile.good()) {
1821     throw DataException("archiveData Error: problem writing data to archive file");
1822     }
1823    
1824     //
1825     // Close archive file
1826     archiveFile.close();
1827    
1828     if (!archiveFile.good()) {
1829     throw DataException("archiveData Error: problem closing archive file");
1830     }
1831    
1832 jgs 119 }
1833    
1834     void
1835     Data::extractData(const std::string fileName,
1836     const FunctionSpace& fspace)
1837     {
1838     //
1839     // Can only extract Data to an object which is initially DataEmpty
1840     if (!isEmpty()) {
1841     throw DataException("extractData Error: can only extract to DataEmpty object");
1842     }
1843    
1844     cout << "Extracting Data object from: " << fileName << endl;
1845    
1846     int dataType;
1847     int noSamples;
1848     int noDPPSample;
1849     int functionSpaceType;
1850     int dataPointRank;
1851     int dataPointSize;
1852     int dataLength;
1853     DataArrayView::ShapeType dataPointShape;
1854     int flatShape[4];
1855    
1856     //
1857 jgs 123 // Open the archive file
1858 jgs 119 ifstream archiveFile;
1859     archiveFile.open(fileName.data(), ios::in);
1860    
1861     if (!archiveFile.good()) {
1862     throw DataException("extractData Error: problem opening archive file");
1863     }
1864    
1865 jgs 123 //
1866     // Read common data items from archive file
1867 jgs 119 archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1868     archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1869     archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1870     archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1871     archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1872     archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1873     archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1874     for (int dim = 0; dim < 4; dim++) {
1875     archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1876     if (flatShape[dim]>0) {
1877     dataPointShape.push_back(flatShape[dim]);
1878     }
1879     }
1880     int referenceNumbers[noSamples];
1881     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1882     archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1883     }
1884     int tagNumbers[noSamples];
1885     if (dataType==2) {
1886     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1887     archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1888     }
1889     }
1890    
1891     if (!archiveFile.good()) {
1892     throw DataException("extractData Error: problem reading from archive file");
1893     }
1894    
1895 jgs 123 //
1896     // Verify the values just read from the archive file
1897 jgs 119 switch (dataType) {
1898     case 0:
1899     cout << "\tdataType: DataEmpty" << endl;
1900     break;
1901     case 1:
1902     cout << "\tdataType: DataConstant" << endl;
1903     break;
1904     case 2:
1905     cout << "\tdataType: DataTagged" << endl;
1906     break;
1907     case 3:
1908     cout << "\tdataType: DataExpanded" << endl;
1909     break;
1910     default:
1911     throw DataException("extractData Error: undefined dataType read from archive file");
1912     break;
1913     }
1914    
1915     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1916     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1917     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1918     cout << "\tshape: < ";
1919     for (int dim = 0; dim < dataPointRank; dim++) {
1920     cout << dataPointShape[dim] << " ";
1921     }
1922     cout << ">" << endl;
1923    
1924     //
1925     // Verify that supplied FunctionSpace object is compatible with this Data object.
1926     if ( (fspace.getTypeCode()!=functionSpaceType) ||
1927     (fspace.getNumSamples()!=noSamples) ||
1928     (fspace.getNumDPPSample()!=noDPPSample)
1929     ) {
1930     throw DataException("extractData Error: incompatible FunctionSpace");
1931     }
1932     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1933     if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
1934     throw DataException("extractData Error: incompatible FunctionSpace");
1935     }
1936     }
1937     if (dataType==2) {
1938     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1939     if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
1940     throw DataException("extractData Error: incompatible FunctionSpace");
1941     }
1942     }
1943     }
1944    
1945     //
1946     // Construct a DataVector to hold underlying data values
1947     DataVector dataVec(dataLength);
1948    
1949     //
1950     // Load this DataVector with the appropriate values
1951 jgs 123 int noValues;
1952     archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
1953     cout << "\tnoValues: " << noValues << endl;
1954 jgs 119 switch (dataType) {
1955     case 0:
1956     // DataEmpty
1957 jgs 123 if (noValues != 0) {
1958     throw DataException("extractData Error: problem reading data from archive file");
1959     }
1960 jgs 119 break;
1961     case 1:
1962     // DataConstant
1963 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
1964     throw DataException("extractData Error: problem reading data from archive file");
1965     }
1966 jgs 119 break;
1967     case 2:
1968     // DataTagged
1969 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
1970     throw DataException("extractData Error: problem reading data from archive file");
1971     }
1972 jgs 119 break;
1973     case 3:
1974     // DataExpanded
1975 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
1976     throw DataException("extractData Error: problem reading data from archive file");
1977     }
1978 jgs 119 break;
1979     }
1980    
1981 jgs 123 if (!archiveFile.good()) {
1982     throw DataException("extractData Error: problem reading from archive file");
1983     }
1984    
1985 jgs 119 //
1986 jgs 123 // Close archive file
1987     archiveFile.close();
1988    
1989     if (!archiveFile.good()) {
1990     throw DataException("extractData Error: problem closing archive file");
1991     }
1992    
1993     //
1994 jgs 119 // Construct an appropriate Data object
1995     DataAbstract* tempData;
1996     switch (dataType) {
1997     case 0:
1998     // DataEmpty
1999     tempData=new DataEmpty();
2000     break;
2001     case 1:
2002     // DataConstant
2003     tempData=new DataConstant(fspace,dataPointShape,dataVec);
2004     break;
2005     case 2:
2006     // DataTagged
2007     tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2008     break;
2009     case 3:
2010     // DataExpanded
2011     tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2012     break;
2013     }
2014     shared_ptr<DataAbstract> temp_data(tempData);
2015     m_data=temp_data;
2016     }
2017    
2018 jgs 94 ostream& escript::operator<<(ostream& o, const Data& data)
2019     {
2020     o << data.toString();
2021     return o;
2022     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26