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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 126 - (hide annotations)
Fri Jul 22 03:53:08 2005 UTC (13 years, 10 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/Data.cpp
File size: 51343 byte(s)
Merge of development branch back to main trunk on 2005-07-22

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26