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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 775 - (hide annotations)
Mon Jul 10 04:00:08 2006 UTC (13 years, 4 months ago) by ksteube
File size: 59356 byte(s)
Modified the following python methods in escript/py_src/util.py to
call faster C++ methods:
	escript_trace
	escript_transpose
	escript_symmetric
	escript_nonsymmetric

1 jgs 94 // $Id$
2 jgs 480
3 elspeth 615 /*
4     ************************************************************
5     * Copyright 2006 by ACcESS MNRF *
6     * *
7     * http://www.access.edu.au *
8     * Primary Business: Queensland, Australia *
9     * Licensed under the Open Software License version 3.0 *
10     * http://www.opensource.org/licenses/osl-3.0.php *
11     * *
12     ************************************************************
13     */
14 jgs 474 #include "Data.h"
15 jgs 94
16 jgs 480 #include "DataExpanded.h"
17     #include "DataConstant.h"
18     #include "DataTagged.h"
19     #include "DataEmpty.h"
20     #include "DataArray.h"
21     #include "DataArrayView.h"
22     #include "DataProf.h"
23     #include "FunctionSpaceFactory.h"
24     #include "AbstractContinuousDomain.h"
25     #include "UnaryFuncs.h"
26    
27 jgs 119 #include <fstream>
28 jgs 94 #include <algorithm>
29     #include <vector>
30     #include <functional>
31    
32 jgs 153 #include <boost/python/dict.hpp>
33 jgs 94 #include <boost/python/extract.hpp>
34     #include <boost/python/long.hpp>
35    
36     using namespace std;
37     using namespace boost::python;
38     using namespace boost;
39     using namespace escript;
40    
41 jgs 151 #if defined DOPROF
42 jgs 123 //
43     // global table of profiling data for all Data objects
44     DataProf dataProfTable;
45 jgs 151 #endif
46 jgs 123
47 jgs 94 Data::Data()
48     {
49     //
50     // Default data is type DataEmpty
51     DataAbstract* temp=new DataEmpty();
52 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
53     m_data=temp_data;
54 jgs 151 #if defined DOPROF
55 jgs 123 // create entry in global profiling table for this object
56     profData = dataProfTable.newData();
57 jgs 151 #endif
58 jgs 94 }
59    
60     Data::Data(double value,
61     const tuple& shape,
62     const FunctionSpace& what,
63     bool expanded)
64     {
65     DataArrayView::ShapeType dataPointShape;
66     for (int i = 0; i < shape.attr("__len__")(); ++i) {
67     dataPointShape.push_back(extract<const int>(shape[i]));
68     }
69     DataArray temp(dataPointShape,value);
70     initialise(temp.getView(),what,expanded);
71 jgs 151 #if defined DOPROF
72 jgs 123 // create entry in global profiling table for this object
73     profData = dataProfTable.newData();
74 jgs 151 #endif
75 jgs 94 }
76    
77     Data::Data(double value,
78     const DataArrayView::ShapeType& dataPointShape,
79     const FunctionSpace& what,
80     bool expanded)
81     {
82     DataArray temp(dataPointShape,value);
83     pair<int,int> dataShape=what.getDataShape();
84     initialise(temp.getView(),what,expanded);
85 jgs 151 #if defined DOPROF
86 jgs 123 // create entry in global profiling table for this object
87     profData = dataProfTable.newData();
88 jgs 151 #endif
89 jgs 94 }
90    
91 jgs 102 Data::Data(const Data& inData)
92 jgs 94 {
93 jgs 102 m_data=inData.m_data;
94 jgs 151 #if defined DOPROF
95 jgs 123 // create entry in global profiling table for this object
96     profData = dataProfTable.newData();
97 jgs 151 #endif
98 jgs 94 }
99    
100     Data::Data(const Data& inData,
101     const DataArrayView::RegionType& region)
102     {
103     //
104 jgs 102 // Create Data which is a slice of another Data
105     DataAbstract* tmp = inData.m_data->getSlice(region);
106     shared_ptr<DataAbstract> temp_data(tmp);
107     m_data=temp_data;
108 jgs 151 #if defined DOPROF
109 jgs 123 // create entry in global profiling table for this object
110     profData = dataProfTable.newData();
111 jgs 151 #endif
112 jgs 94 }
113    
114     Data::Data(const Data& inData,
115     const FunctionSpace& functionspace)
116     {
117 gross 711 #if defined DOPROF
118     // create entry in global profiling table for this object
119     profData = dataProfTable.newData();
120     #endif
121 jgs 94 if (inData.getFunctionSpace()==functionspace) {
122     m_data=inData.m_data;
123     } else {
124 gross 711 #if defined DOPROF
125     profData->interpolate++;
126     #endif
127 jgs 94 Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
128 jgs 123 // Note: Must use a reference or pointer to a derived object
129 jgs 94 // in order to get polymorphic behaviour. Shouldn't really
130     // be able to create an instance of AbstractDomain but that was done
131 jgs 123 // as a boost:python work around which may no longer be required.
132 jgs 94 const AbstractDomain& inDataDomain=inData.getDomain();
133     if (inDataDomain==functionspace.getDomain()) {
134     inDataDomain.interpolateOnDomain(tmp,inData);
135     } else {
136     inDataDomain.interpolateACross(tmp,inData);
137     }
138     m_data=tmp.m_data;
139     }
140     }
141    
142     Data::Data(const DataTagged::TagListType& tagKeys,
143     const DataTagged::ValueListType & values,
144     const DataArrayView& defaultValue,
145     const FunctionSpace& what,
146     bool expanded)
147     {
148     DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
149 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
150     m_data=temp_data;
151 jgs 94 if (expanded) {
152     expand();
153     }
154 jgs 151 #if defined DOPROF
155 jgs 123 // create entry in global profiling table for this object
156     profData = dataProfTable.newData();
157 jgs 151 #endif
158 jgs 94 }
159    
160     Data::Data(const numeric::array& value,
161     const FunctionSpace& what,
162     bool expanded)
163     {
164     initialise(value,what,expanded);
165 jgs 151 #if defined DOPROF
166 jgs 123 // create entry in global profiling table for this object
167     profData = dataProfTable.newData();
168 jgs 151 #endif
169 jgs 94 }
170    
171     Data::Data(const DataArrayView& value,
172     const FunctionSpace& what,
173     bool expanded)
174     {
175     initialise(value,what,expanded);
176 jgs 151 #if defined DOPROF
177 jgs 123 // create entry in global profiling table for this object
178     profData = dataProfTable.newData();
179 jgs 151 #endif
180 jgs 94 }
181    
182     Data::Data(const object& value,
183     const FunctionSpace& what,
184     bool expanded)
185     {
186     numeric::array asNumArray(value);
187     initialise(asNumArray,what,expanded);
188 jgs 151 #if defined DOPROF
189 jgs 123 // create entry in global profiling table for this object
190     profData = dataProfTable.newData();
191 jgs 151 #endif
192 jgs 94 }
193    
194     Data::Data(const object& value,
195     const Data& other)
196     {
197     //
198     // Create DataConstant using the given value and all other parameters
199     // copied from other. If value is a rank 0 object this Data
200     // will assume the point data shape of other.
201     DataArray temp(value);
202     if (temp.getView().getRank()==0) {
203     //
204     // Create a DataArray with the scalar value for all elements
205     DataArray temp2(other.getPointDataView().getShape(),temp.getView()());
206     initialise(temp2.getView(),other.getFunctionSpace(),false);
207     } else {
208     //
209     // Create a DataConstant with the same sample shape as other
210     initialise(temp.getView(),other.getFunctionSpace(),false);
211     }
212 jgs 151 #if defined DOPROF
213 jgs 123 // create entry in global profiling table for this object
214     profData = dataProfTable.newData();
215 jgs 151 #endif
216 jgs 94 }
217    
218 jgs 151 Data::~Data()
219     {
220    
221     }
222    
223 jgs 94 escriptDataC
224     Data::getDataC()
225     {
226     escriptDataC temp;
227     temp.m_dataPtr=(void*)this;
228     return temp;
229     }
230    
231     escriptDataC
232     Data::getDataC() const
233     {
234     escriptDataC temp;
235     temp.m_dataPtr=(void*)this;
236     return temp;
237     }
238    
239 jgs 121 const boost::python::tuple
240 jgs 94 Data::getShapeTuple() const
241     {
242     const DataArrayView::ShapeType& shape=getDataPointShape();
243     switch(getDataPointRank()) {
244     case 0:
245     return make_tuple();
246     case 1:
247     return make_tuple(long_(shape[0]));
248     case 2:
249     return make_tuple(long_(shape[0]),long_(shape[1]));
250     case 3:
251     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]));
252     case 4:
253     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]),long_(shape[3]));
254     default:
255     throw DataException("Error - illegal Data rank.");
256     }
257     }
258    
259     void
260     Data::copy(const Data& other)
261     {
262     //
263     // Perform a deep copy
264     {
265     DataExpanded* temp=dynamic_cast<DataExpanded*>(other.m_data.get());
266     if (temp!=0) {
267     //
268     // Construct a DataExpanded copy
269     DataAbstract* newData=new DataExpanded(*temp);
270 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
271     m_data=temp_data;
272 jgs 94 return;
273     }
274     }
275     {
276     DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());
277     if (temp!=0) {
278     //
279 jgs 102 // Construct a DataTagged copy
280 jgs 94 DataAbstract* newData=new DataTagged(*temp);
281 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
282     m_data=temp_data;
283 jgs 94 return;
284     }
285     }
286     {
287     DataConstant* temp=dynamic_cast<DataConstant*>(other.m_data.get());
288     if (temp!=0) {
289     //
290     // Construct a DataConstant copy
291     DataAbstract* newData=new DataConstant(*temp);
292 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
293     m_data=temp_data;
294 jgs 94 return;
295     }
296     }
297 jgs 102 {
298     DataEmpty* temp=dynamic_cast<DataEmpty*>(other.m_data.get());
299     if (temp!=0) {
300     //
301     // Construct a DataEmpty copy
302     DataAbstract* newData=new DataEmpty();
303     shared_ptr<DataAbstract> temp_data(newData);
304     m_data=temp_data;
305     return;
306     }
307     }
308 jgs 94 throw DataException("Error - Copy not implemented for this Data type.");
309     }
310    
311     void
312     Data::copyWithMask(const Data& other,
313     const Data& mask)
314     {
315     Data mask1;
316     Data mask2;
317    
318     mask1 = mask.wherePositive();
319     mask2.copy(mask1);
320    
321     mask1 *= other;
322     mask2 *= *this;
323     mask2 = *this - mask2;
324    
325     *this = mask1 + mask2;
326     }
327    
328     bool
329     Data::isExpanded() const
330     {
331     DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
332     return (temp!=0);
333     }
334    
335     bool
336     Data::isTagged() const
337     {
338     DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
339     return (temp!=0);
340     }
341    
342 bcumming 751 /* TODO */
343     /* global reduction -- the local data being empty does not imply that it is empty on other processers*/
344 jgs 94 bool
345     Data::isEmpty() const
346     {
347     DataEmpty* temp=dynamic_cast<DataEmpty*>(m_data.get());
348     return (temp!=0);
349     }
350    
351     bool
352     Data::isConstant() const
353     {
354     DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
355     return (temp!=0);
356     }
357    
358     void
359     Data::expand()
360     {
361     if (isConstant()) {
362     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
363     DataAbstract* temp=new DataExpanded(*tempDataConst);
364 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
365     m_data=temp_data;
366 jgs 94 } else if (isTagged()) {
367     DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
368     DataAbstract* temp=new DataExpanded(*tempDataTag);
369 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
370     m_data=temp_data;
371 jgs 94 } else if (isExpanded()) {
372     //
373     // do nothing
374     } else if (isEmpty()) {
375     throw DataException("Error - Expansion of DataEmpty not possible.");
376     } else {
377     throw DataException("Error - Expansion not implemented for this Data type.");
378     }
379     }
380    
381     void
382     Data::tag()
383     {
384     if (isConstant()) {
385     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
386     DataAbstract* temp=new DataTagged(*tempDataConst);
387 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
388     m_data=temp_data;
389 jgs 94 } else if (isTagged()) {
390     // do nothing
391     } else if (isExpanded()) {
392     throw DataException("Error - Creating tag data from DataExpanded not possible.");
393     } else if (isEmpty()) {
394     throw DataException("Error - Creating tag data from DataEmpty not possible.");
395     } else {
396     throw DataException("Error - Tagging not implemented for this Data type.");
397     }
398     }
399    
400 jgs 102 void
401     Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)
402     {
403     m_data->reshapeDataPoint(shape);
404     }
405    
406 jgs 94 Data
407 gross 698 Data::wherePositive() const
408 jgs 94 {
409 jgs 151 #if defined DOPROF
410 jgs 123 profData->where++;
411 jgs 151 #endif
412 gross 698 return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
413 jgs 94 }
414    
415     Data
416 gross 698 Data::whereNegative() const
417 jgs 102 {
418 jgs 151 #if defined DOPROF
419 jgs 123 profData->where++;
420 jgs 151 #endif
421 gross 698 return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
422 jgs 102 }
423    
424     Data
425 gross 698 Data::whereNonNegative() const
426 jgs 94 {
427 jgs 151 #if defined DOPROF
428 jgs 123 profData->where++;
429 jgs 151 #endif
430 gross 698 return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
431 jgs 94 }
432    
433     Data
434 gross 698 Data::whereNonPositive() const
435 jgs 94 {
436 jgs 151 #if defined DOPROF
437 jgs 123 profData->where++;
438 jgs 151 #endif
439 gross 698 return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
440 jgs 94 }
441    
442     Data
443 jgs 571 Data::whereZero(double tol) const
444 jgs 94 {
445 jgs 151 #if defined DOPROF
446 jgs 123 profData->where++;
447 jgs 151 #endif
448 jgs 571 Data dataAbs=abs();
449     return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));
450 jgs 94 }
451    
452     Data
453 jgs 571 Data::whereNonZero(double tol) const
454 jgs 102 {
455 jgs 151 #if defined DOPROF
456 jgs 123 profData->where++;
457 jgs 151 #endif
458 jgs 571 Data dataAbs=abs();
459     return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));
460 jgs 102 }
461    
462     Data
463 jgs 94 Data::interpolate(const FunctionSpace& functionspace) const
464     {
465 jgs 151 #if defined DOPROF
466 jgs 123 profData->interpolate++;
467 jgs 151 #endif
468 jgs 94 return Data(*this,functionspace);
469     }
470    
471     bool
472     Data::probeInterpolation(const FunctionSpace& functionspace) const
473     {
474     if (getFunctionSpace()==functionspace) {
475     return true;
476     } else {
477     const AbstractDomain& domain=getDomain();
478     if (domain==functionspace.getDomain()) {
479     return domain.probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());
480     } else {
481     return domain.probeInterpolationACross(getFunctionSpace().getTypeCode(),functionspace.getDomain(),functionspace.getTypeCode());
482     }
483     }
484     }
485    
486     Data
487     Data::gradOn(const FunctionSpace& functionspace) const
488     {
489 jgs 151 #if defined DOPROF
490 jgs 123 profData->grad++;
491 jgs 151 #endif
492 jgs 94 if (functionspace.getDomain()!=getDomain())
493     throw DataException("Error - gradient cannot be calculated on different domains.");
494     DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
495     grad_shape.push_back(functionspace.getDim());
496     Data out(0.0,grad_shape,functionspace,true);
497     getDomain().setToGradient(out,*this);
498     return out;
499     }
500    
501     Data
502     Data::grad() const
503     {
504     return gradOn(escript::function(getDomain()));
505     }
506    
507     int
508     Data::getDataPointSize() const
509     {
510     return getPointDataView().noValues();
511     }
512    
513     DataArrayView::ValueType::size_type
514     Data::getLength() const
515     {
516     return m_data->getLength();
517     }
518    
519     const DataArrayView::ShapeType&
520     Data::getDataPointShape() const
521     {
522     return getPointDataView().getShape();
523     }
524    
525 jgs 126 void
526     Data::fillFromNumArray(const boost::python::numeric::array num_array)
527     {
528     //
529 jgs 149 // check rank
530 jgs 126 if (num_array.getrank()<getDataPointRank())
531     throw DataException("Rank of numarray does not match Data object rank");
532 jgs 149
533 jgs 126 //
534 jgs 149 // check shape of num_array
535 jgs 126 for (int i=0; i<getDataPointRank(); i++) {
536     if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])
537     throw DataException("Shape of numarray does not match Data object rank");
538     }
539 jgs 149
540 jgs 126 //
541     // make sure data is expanded:
542 jgs 149 if (!isExpanded()) {
543     expand();
544     }
545    
546 jgs 126 //
547 jgs 149 // and copy over
548 jgs 126 m_data->copyAll(num_array);
549     }
550    
551 jgs 121 const
552 jgs 94 boost::python::numeric::array
553 jgs 117 Data::convertToNumArray()
554     {
555     //
556 jgs 121 // determine the total number of data points
557 jgs 117 int numSamples = getNumSamples();
558     int numDataPointsPerSample = getNumDataPointsPerSample();
559     int numDataPoints = numSamples * numDataPointsPerSample;
560    
561     //
562     // determine the rank and shape of each data point
563     int dataPointRank = getDataPointRank();
564     DataArrayView::ShapeType dataPointShape = getDataPointShape();
565    
566     //
567     // create the numeric array to be returned
568     boost::python::numeric::array numArray(0.0);
569    
570     //
571     // the rank of the returned numeric array will be the rank of
572     // the data points, plus one. Where the rank of the array is n,
573     // the last n-1 dimensions will be equal to the shape of the
574     // data points, whilst the first dimension will be equal to the
575     // total number of data points. Thus the array will consist of
576     // a serial vector of the data points.
577     int arrayRank = dataPointRank + 1;
578     DataArrayView::ShapeType arrayShape;
579     arrayShape.push_back(numDataPoints);
580     for (int d=0; d<dataPointRank; d++) {
581     arrayShape.push_back(dataPointShape[d]);
582     }
583    
584     //
585     // resize the numeric array to the shape just calculated
586     if (arrayRank==1) {
587     numArray.resize(arrayShape[0]);
588     }
589     if (arrayRank==2) {
590     numArray.resize(arrayShape[0],arrayShape[1]);
591     }
592     if (arrayRank==3) {
593     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
594     }
595     if (arrayRank==4) {
596     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
597     }
598     if (arrayRank==5) {
599     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
600     }
601    
602     //
603     // loop through each data point in turn, loading the values for that data point
604     // into the numeric array.
605     int dataPoint = 0;
606     for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
607     for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
608     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
609     if (dataPointRank==0) {
610     numArray[dataPoint]=dataPointView();
611     }
612     if (dataPointRank==1) {
613     for (int i=0; i<dataPointShape[0]; i++) {
614     numArray[dataPoint][i]=dataPointView(i);
615     }
616     }
617     if (dataPointRank==2) {
618     for (int i=0; i<dataPointShape[0]; i++) {
619     for (int j=0; j<dataPointShape[1]; j++) {
620     numArray[dataPoint][i][j] = dataPointView(i,j);
621     }
622     }
623     }
624     if (dataPointRank==3) {
625     for (int i=0; i<dataPointShape[0]; i++) {
626     for (int j=0; j<dataPointShape[1]; j++) {
627     for (int k=0; k<dataPointShape[2]; k++) {
628     numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
629     }
630     }
631     }
632     }
633     if (dataPointRank==4) {
634     for (int i=0; i<dataPointShape[0]; i++) {
635     for (int j=0; j<dataPointShape[1]; j++) {
636     for (int k=0; k<dataPointShape[2]; k++) {
637     for (int l=0; l<dataPointShape[3]; l++) {
638     numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
639     }
640     }
641     }
642     }
643     }
644     dataPoint++;
645 jgs 121 }
646     }
647 jgs 117
648 jgs 121 //
649     // return the loaded array
650     return numArray;
651     }
652    
653     const
654     boost::python::numeric::array
655     Data::convertToNumArrayFromSampleNo(int sampleNo)
656     {
657     //
658     // Check a valid sample number has been supplied
659     if (sampleNo >= getNumSamples()) {
660     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
661     }
662    
663     //
664     // determine the number of data points per sample
665     int numDataPointsPerSample = getNumDataPointsPerSample();
666    
667     //
668     // determine the rank and shape of each data point
669     int dataPointRank = getDataPointRank();
670     DataArrayView::ShapeType dataPointShape = getDataPointShape();
671    
672     //
673     // create the numeric array to be returned
674     boost::python::numeric::array numArray(0.0);
675    
676     //
677     // the rank of the returned numeric array will be the rank of
678     // the data points, plus one. Where the rank of the array is n,
679     // the last n-1 dimensions will be equal to the shape of the
680     // data points, whilst the first dimension will be equal to the
681     // total number of data points. Thus the array will consist of
682     // a serial vector of the data points.
683     int arrayRank = dataPointRank + 1;
684     DataArrayView::ShapeType arrayShape;
685     arrayShape.push_back(numDataPointsPerSample);
686     for (int d=0; d<dataPointRank; d++) {
687     arrayShape.push_back(dataPointShape[d]);
688     }
689    
690     //
691     // resize the numeric array to the shape just calculated
692     if (arrayRank==1) {
693     numArray.resize(arrayShape[0]);
694     }
695     if (arrayRank==2) {
696     numArray.resize(arrayShape[0],arrayShape[1]);
697     }
698     if (arrayRank==3) {
699     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
700     }
701     if (arrayRank==4) {
702     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
703     }
704     if (arrayRank==5) {
705     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
706     }
707    
708     //
709     // loop through each data point in turn, loading the values for that data point
710     // into the numeric array.
711     for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {
712     DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);
713     if (dataPointRank==0) {
714     numArray[dataPoint]=dataPointView();
715 jgs 117 }
716 jgs 121 if (dataPointRank==1) {
717     for (int i=0; i<dataPointShape[0]; i++) {
718     numArray[dataPoint][i]=dataPointView(i);
719     }
720     }
721     if (dataPointRank==2) {
722     for (int i=0; i<dataPointShape[0]; i++) {
723     for (int j=0; j<dataPointShape[1]; j++) {
724     numArray[dataPoint][i][j] = dataPointView(i,j);
725     }
726     }
727     }
728     if (dataPointRank==3) {
729     for (int i=0; i<dataPointShape[0]; i++) {
730     for (int j=0; j<dataPointShape[1]; j++) {
731     for (int k=0; k<dataPointShape[2]; k++) {
732     numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
733     }
734     }
735     }
736     }
737     if (dataPointRank==4) {
738     for (int i=0; i<dataPointShape[0]; i++) {
739     for (int j=0; j<dataPointShape[1]; j++) {
740     for (int k=0; k<dataPointShape[2]; k++) {
741     for (int l=0; l<dataPointShape[3]; l++) {
742     numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
743     }
744     }
745     }
746     }
747     }
748 jgs 117 }
749    
750     //
751     // return the loaded array
752     return numArray;
753     }
754    
755 jgs 121 const
756 jgs 117 boost::python::numeric::array
757 jgs 121 Data::convertToNumArrayFromDPNo(int sampleNo,
758     int dataPointNo)
759     {
760     //
761     // Check a valid sample number has been supplied
762     if (sampleNo >= getNumSamples()) {
763     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
764     }
765    
766     //
767     // Check a valid data point number has been supplied
768     if (dataPointNo >= getNumDataPointsPerSample()) {
769     throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");
770     }
771    
772     //
773     // determine the rank and shape of each data point
774     int dataPointRank = getDataPointRank();
775     DataArrayView::ShapeType dataPointShape = getDataPointShape();
776    
777     //
778     // create the numeric array to be returned
779     boost::python::numeric::array numArray(0.0);
780    
781     //
782     // the shape of the returned numeric array will be the same
783     // as that of the data point
784     int arrayRank = dataPointRank;
785     DataArrayView::ShapeType arrayShape = dataPointShape;
786    
787     //
788     // resize the numeric array to the shape just calculated
789     if (arrayRank==0) {
790     numArray.resize(1);
791     }
792     if (arrayRank==1) {
793     numArray.resize(arrayShape[0]);
794     }
795     if (arrayRank==2) {
796     numArray.resize(arrayShape[0],arrayShape[1]);
797     }
798     if (arrayRank==3) {
799     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
800     }
801     if (arrayRank==4) {
802     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
803     }
804    
805     //
806     // load the values for the data point into the numeric array.
807     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
808     if (dataPointRank==0) {
809     numArray[0]=dataPointView();
810     }
811     if (dataPointRank==1) {
812     for (int i=0; i<dataPointShape[0]; i++) {
813     numArray[i]=dataPointView(i);
814     }
815     }
816     if (dataPointRank==2) {
817     for (int i=0; i<dataPointShape[0]; i++) {
818     for (int j=0; j<dataPointShape[1]; j++) {
819     numArray[i][j] = dataPointView(i,j);
820     }
821     }
822     }
823     if (dataPointRank==3) {
824     for (int i=0; i<dataPointShape[0]; i++) {
825     for (int j=0; j<dataPointShape[1]; j++) {
826     for (int k=0; k<dataPointShape[2]; k++) {
827     numArray[i][j][k]=dataPointView(i,j,k);
828     }
829     }
830     }
831     }
832     if (dataPointRank==4) {
833     for (int i=0; i<dataPointShape[0]; i++) {
834     for (int j=0; j<dataPointShape[1]; j++) {
835     for (int k=0; k<dataPointShape[2]; k++) {
836     for (int l=0; l<dataPointShape[3]; l++) {
837     numArray[i][j][k][l]=dataPointView(i,j,k,l);
838     }
839     }
840     }
841     }
842     }
843    
844     //
845     // return the loaded array
846     return numArray;
847     }
848    
849     boost::python::numeric::array
850 jgs 94 Data::integrate() const
851     {
852     int index;
853     int rank = getDataPointRank();
854     DataArrayView::ShapeType shape = getDataPointShape();
855    
856 jgs 151 #if defined DOPROF
857 jgs 123 profData->integrate++;
858 jgs 151 #endif
859 jgs 123
860 jgs 94 //
861     // calculate the integral values
862     vector<double> integrals(getDataPointSize());
863     AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);
864    
865     //
866     // create the numeric array to be returned
867     // and load the array with the integral values
868     boost::python::numeric::array bp_array(1.0);
869     if (rank==0) {
870 jgs 108 bp_array.resize(1);
871 jgs 94 index = 0;
872     bp_array[0] = integrals[index];
873     }
874     if (rank==1) {
875     bp_array.resize(shape[0]);
876     for (int i=0; i<shape[0]; i++) {
877     index = i;
878     bp_array[i] = integrals[index];
879     }
880     }
881     if (rank==2) {
882 gross 436 bp_array.resize(shape[0],shape[1]);
883     for (int i=0; i<shape[0]; i++) {
884     for (int j=0; j<shape[1]; j++) {
885     index = i + shape[0] * j;
886     bp_array[make_tuple(i,j)] = integrals[index];
887     }
888     }
889 jgs 94 }
890     if (rank==3) {
891     bp_array.resize(shape[0],shape[1],shape[2]);
892     for (int i=0; i<shape[0]; i++) {
893     for (int j=0; j<shape[1]; j++) {
894     for (int k=0; k<shape[2]; k++) {
895     index = i + shape[0] * ( j + shape[1] * k );
896 gross 436 bp_array[make_tuple(i,j,k)] = integrals[index];
897 jgs 94 }
898     }
899     }
900     }
901     if (rank==4) {
902     bp_array.resize(shape[0],shape[1],shape[2],shape[3]);
903     for (int i=0; i<shape[0]; i++) {
904     for (int j=0; j<shape[1]; j++) {
905     for (int k=0; k<shape[2]; k++) {
906     for (int l=0; l<shape[3]; l++) {
907     index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
908 gross 436 bp_array[make_tuple(i,j,k,l)] = integrals[index];
909 jgs 94 }
910     }
911     }
912     }
913     }
914    
915     //
916     // return the loaded array
917     return bp_array;
918     }
919    
920     Data
921     Data::sin() const
922     {
923 jgs 151 #if defined DOPROF
924 jgs 123 profData->unary++;
925 jgs 151 #endif
926 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
927     }
928    
929     Data
930     Data::cos() const
931     {
932 jgs 151 #if defined DOPROF
933 jgs 123 profData->unary++;
934 jgs 151 #endif
935 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
936     }
937    
938     Data
939     Data::tan() const
940     {
941 jgs 151 #if defined DOPROF
942 jgs 123 profData->unary++;
943 jgs 151 #endif
944 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
945     }
946    
947     Data
948 jgs 150 Data::asin() const
949     {
950 jgs 151 #if defined DOPROF
951 jgs 150 profData->unary++;
952 jgs 151 #endif
953 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
954     }
955    
956     Data
957     Data::acos() const
958     {
959 jgs 151 #if defined DOPROF
960 jgs 150 profData->unary++;
961 jgs 151 #endif
962 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
963     }
964    
965     Data
966     Data::atan() const
967     {
968 jgs 151 #if defined DOPROF
969 jgs 150 profData->unary++;
970 jgs 151 #endif
971 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
972     }
973    
974     Data
975     Data::sinh() const
976     {
977 jgs 151 #if defined DOPROF
978 jgs 150 profData->unary++;
979 jgs 151 #endif
980 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
981     }
982    
983     Data
984     Data::cosh() const
985     {
986 jgs 151 #if defined DOPROF
987 jgs 150 profData->unary++;
988 jgs 151 #endif
989 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
990     }
991    
992     Data
993     Data::tanh() const
994     {
995 jgs 151 #if defined DOPROF
996 jgs 150 profData->unary++;
997 jgs 151 #endif
998 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
999     }
1000    
1001     Data
1002     Data::asinh() const
1003     {
1004 jgs 151 #if defined DOPROF
1005 jgs 150 profData->unary++;
1006 jgs 151 #endif
1007 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
1008     }
1009    
1010     Data
1011     Data::acosh() const
1012     {
1013 jgs 151 #if defined DOPROF
1014 jgs 150 profData->unary++;
1015 jgs 151 #endif
1016 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
1017     }
1018    
1019     Data
1020     Data::atanh() const
1021     {
1022 jgs 151 #if defined DOPROF
1023 jgs 150 profData->unary++;
1024 jgs 151 #endif
1025 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
1026     }
1027    
1028     Data
1029 gross 286 Data::log10() const
1030 jgs 94 {
1031 jgs 151 #if defined DOPROF
1032 jgs 123 profData->unary++;
1033 jgs 151 #endif
1034 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
1035     }
1036    
1037     Data
1038 gross 286 Data::log() const
1039 jgs 94 {
1040 jgs 151 #if defined DOPROF
1041 jgs 123 profData->unary++;
1042 jgs 151 #endif
1043 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
1044     }
1045    
1046 jgs 106 Data
1047     Data::sign() const
1048 jgs 94 {
1049 jgs 151 #if defined DOPROF
1050 jgs 123 profData->unary++;
1051 jgs 151 #endif
1052 jgs 106 return escript::unaryOp(*this,escript::fsign);
1053 jgs 94 }
1054    
1055 jgs 106 Data
1056     Data::abs() const
1057 jgs 94 {
1058 jgs 151 #if defined DOPROF
1059 jgs 123 profData->unary++;
1060 jgs 151 #endif
1061 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
1062 jgs 94 }
1063    
1064 jgs 106 Data
1065     Data::neg() const
1066 jgs 94 {
1067 jgs 151 #if defined DOPROF
1068 jgs 123 profData->unary++;
1069 jgs 151 #endif
1070 jgs 106 return escript::unaryOp(*this,negate<double>());
1071 jgs 94 }
1072    
1073 jgs 102 Data
1074 jgs 106 Data::pos() const
1075 jgs 94 {
1076 jgs 151 #if defined DOPROF
1077 jgs 123 profData->unary++;
1078 jgs 151 #endif
1079 jgs 148 Data result;
1080     // perform a deep copy
1081     result.copy(*this);
1082     return result;
1083 jgs 102 }
1084    
1085     Data
1086 jgs 106 Data::exp() const
1087 jgs 102 {
1088 jgs 151 #if defined DOPROF
1089 jgs 123 profData->unary++;
1090 jgs 151 #endif
1091 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1092 jgs 102 }
1093    
1094     Data
1095 jgs 106 Data::sqrt() const
1096 jgs 102 {
1097 jgs 151 #if defined DOPROF
1098 jgs 123 profData->unary++;
1099 jgs 151 #endif
1100 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1101 jgs 102 }
1102    
1103 jgs 106 double
1104     Data::Lsup() const
1105 jgs 102 {
1106 bcumming 751 double localValue, globalValue;
1107 jgs 151 #if defined DOPROF
1108 jgs 123 profData->reduction1++;
1109 jgs 151 #endif
1110 jgs 106 //
1111     // set the initial absolute maximum value to zero
1112 bcumming 751
1113 jgs 147 AbsMax abs_max_func;
1114 bcumming 751 localValue = algorithm(abs_max_func,0);
1115     #ifdef PASO_MPI
1116     MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1117     return globalValue;
1118     #else
1119     return localValue;
1120     #endif
1121 jgs 102 }
1122    
1123 jgs 106 double
1124 jgs 117 Data::Linf() const
1125     {
1126 bcumming 751 double localValue, globalValue;
1127 jgs 151 #if defined DOPROF
1128 jgs 123 profData->reduction1++;
1129 jgs 151 #endif
1130 jgs 117 //
1131     // set the initial absolute minimum value to max double
1132 jgs 147 AbsMin abs_min_func;
1133 bcumming 751 localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1134    
1135     #ifdef PASO_MPI
1136     MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1137     return globalValue;
1138     #else
1139     return localValue;
1140     #endif
1141 jgs 117 }
1142    
1143     double
1144 jgs 106 Data::sup() const
1145 jgs 102 {
1146 bcumming 751 double localValue, globalValue;
1147 jgs 151 #if defined DOPROF
1148 jgs 123 profData->reduction1++;
1149 jgs 151 #endif
1150 jgs 106 //
1151     // set the initial maximum value to min possible double
1152 jgs 147 FMax fmax_func;
1153 bcumming 751 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1154     #ifdef PASO_MPI
1155     MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1156     return globalValue;
1157     #else
1158     return localValue;
1159     #endif
1160 jgs 102 }
1161    
1162 jgs 106 double
1163     Data::inf() const
1164 jgs 102 {
1165 bcumming 751 double localValue, globalValue;
1166 jgs 151 #if defined DOPROF
1167 jgs 123 profData->reduction1++;
1168 jgs 151 #endif
1169 jgs 106 //
1170     // set the initial minimum value to max possible double
1171 jgs 147 FMin fmin_func;
1172 bcumming 751 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1173     #ifdef PASO_MPI
1174     MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1175     return globalValue;
1176     #else
1177     return localValue;
1178     #endif
1179 jgs 102 }
1180    
1181 bcumming 751 /* TODO */
1182     /* global reduction */
1183 jgs 102 Data
1184 jgs 106 Data::maxval() const
1185 jgs 102 {
1186 jgs 151 #if defined DOPROF
1187 jgs 123 profData->reduction2++;
1188 jgs 151 #endif
1189 jgs 113 //
1190     // set the initial maximum value to min possible double
1191 jgs 147 FMax fmax_func;
1192     return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1193 jgs 102 }
1194    
1195     Data
1196 jgs 106 Data::minval() const
1197 jgs 102 {
1198 jgs 151 #if defined DOPROF
1199 jgs 123 profData->reduction2++;
1200 jgs 151 #endif
1201 jgs 113 //
1202     // set the initial minimum value to max possible double
1203 jgs 147 FMin fmin_func;
1204     return dp_algorithm(fmin_func,numeric_limits<double>::max());
1205 jgs 102 }
1206    
1207 jgs 123 Data
1208     Data::trace() const
1209     {
1210 jgs 151 #if defined DOPROF
1211 jgs 123 profData->reduction2++;
1212 jgs 151 #endif
1213 jgs 147 Trace trace_func;
1214     return dp_algorithm(trace_func,0);
1215 jgs 123 }
1216    
1217     Data
1218 ksteube 775 Data::symmetric() const
1219 jgs 123 {
1220 ksteube 775 #if defined DOPROF
1221     profData->unary++;
1222     #endif
1223     // check input
1224     DataArrayView::ShapeType s=getDataPointShape();
1225     if (getDataPointRank()==2) {
1226     if(s[0] != s[1])
1227     throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1228     }
1229     else if (getDataPointRank()==4) {
1230     if(!(s[0] == s[2] && s[1] == s[3]))
1231     throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1232     }
1233     else {
1234     throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1235     }
1236     Data ev(0.,getDataPointShape(),getFunctionSpace());
1237     ev.typeMatchRight(*this);
1238     m_data->symmetric(ev.m_data.get());
1239     return ev;
1240     }
1241    
1242     Data
1243     Data::nonsymmetric() const
1244     {
1245     #if defined DOPROF
1246     profData->unary++;
1247     #endif
1248     // check input
1249     DataArrayView::ShapeType s=getDataPointShape();
1250     if (getDataPointRank()==2) {
1251     if(s[0] != s[1])
1252     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1253     DataArrayView::ShapeType ev_shape;
1254     ev_shape.push_back(s[0]);
1255     ev_shape.push_back(s[1]);
1256     Data ev(0.,ev_shape,getFunctionSpace());
1257     ev.typeMatchRight(*this);
1258     m_data->nonsymmetric(ev.m_data.get());
1259     return ev;
1260     }
1261     else if (getDataPointRank()==4) {
1262     if(!(s[0] == s[2] && s[1] == s[3]))
1263     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1264     DataArrayView::ShapeType ev_shape;
1265     ev_shape.push_back(s[0]);
1266     ev_shape.push_back(s[1]);
1267     ev_shape.push_back(s[2]);
1268     ev_shape.push_back(s[3]);
1269     Data ev(0.,ev_shape,getFunctionSpace());
1270     ev.typeMatchRight(*this);
1271     m_data->nonsymmetric(ev.m_data.get());
1272     return ev;
1273     }
1274     else {
1275     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1276     }
1277     }
1278    
1279     Data
1280     Data::matrixtrace(int axis_offset) const
1281     {
1282     #if defined DOPROF
1283     profData->unary++;
1284     #endif
1285     DataArrayView::ShapeType s=getDataPointShape();
1286     if (getDataPointRank()==2) {
1287     DataArrayView::ShapeType ev_shape;
1288     Data ev(0.,ev_shape,getFunctionSpace());
1289     ev.typeMatchRight(*this);
1290     m_data->matrixtrace(ev.m_data.get(), axis_offset);
1291     return ev;
1292     }
1293     if (getDataPointRank()==3) {
1294     DataArrayView::ShapeType ev_shape;
1295     if (axis_offset==0) {
1296     int s2=s[2];
1297     ev_shape.push_back(s2);
1298     }
1299     else if (axis_offset==1) {
1300     int s0=s[0];
1301     ev_shape.push_back(s0);
1302     }
1303     Data ev(0.,ev_shape,getFunctionSpace());
1304     ev.typeMatchRight(*this);
1305     m_data->matrixtrace(ev.m_data.get(), axis_offset);
1306     return ev;
1307     }
1308     if (getDataPointRank()==4) {
1309     DataArrayView::ShapeType ev_shape;
1310     if (axis_offset==0) {
1311     ev_shape.push_back(s[2]);
1312     ev_shape.push_back(s[3]);
1313     }
1314     else if (axis_offset==1) {
1315     ev_shape.push_back(s[0]);
1316     ev_shape.push_back(s[3]);
1317     }
1318     else if (axis_offset==2) {
1319     ev_shape.push_back(s[0]);
1320     ev_shape.push_back(s[1]);
1321     }
1322     Data ev(0.,ev_shape,getFunctionSpace());
1323     ev.typeMatchRight(*this);
1324     m_data->matrixtrace(ev.m_data.get(), axis_offset);
1325     return ev;
1326     }
1327     else {
1328     throw DataException("Error - Data::matrixtrace can only be calculated for rank 2, 3 or 4 object.");
1329     }
1330     }
1331    
1332     Data
1333     Data::transpose(int axis_offset) const
1334     {
1335 jgs 151 #if defined DOPROF
1336 ksteube 775 profData->reduction2++;
1337 jgs 151 #endif
1338 ksteube 775 DataArrayView::ShapeType s=getDataPointShape();
1339     DataArrayView::ShapeType ev_shape;
1340     // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1341     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1342     int rank=getDataPointRank();
1343     if (axis_offset<0 || axis_offset>rank) {
1344     throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1345     }
1346     for (int i=0; i<rank; i++) {
1347     int index = (axis_offset+i)%rank;
1348     ev_shape.push_back(s[index]); // Append to new shape
1349     }
1350     Data ev(0.,ev_shape,getFunctionSpace());
1351     ev.typeMatchRight(*this);
1352     m_data->transpose(ev.m_data.get(), axis_offset);
1353     return ev;
1354 jgs 123 }
1355    
1356 gross 576 Data
1357     Data::eigenvalues() const
1358     {
1359     #if defined DOPROF
1360     profData->unary++;
1361     #endif
1362     // check input
1363     DataArrayView::ShapeType s=getDataPointShape();
1364     if (getDataPointRank()!=2)
1365     throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1366     if(s[0] != s[1])
1367     throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1368     // create return
1369     DataArrayView::ShapeType ev_shape(1,s[0]);
1370     Data ev(0.,ev_shape,getFunctionSpace());
1371     ev.typeMatchRight(*this);
1372     m_data->eigenvalues(ev.m_data.get());
1373     return ev;
1374     }
1375    
1376 jgs 121 const boost::python::tuple
1377 gross 576 Data::eigenvalues_and_eigenvectors(const double tol) const
1378     {
1379     #if defined DOPROF
1380     profData->unary++;
1381     #endif
1382     DataArrayView::ShapeType s=getDataPointShape();
1383     if (getDataPointRank()!=2)
1384     throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1385     if(s[0] != s[1])
1386     throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1387     // create return
1388     DataArrayView::ShapeType ev_shape(1,s[0]);
1389     Data ev(0.,ev_shape,getFunctionSpace());
1390     ev.typeMatchRight(*this);
1391     DataArrayView::ShapeType V_shape(2,s[0]);
1392     Data V(0.,V_shape,getFunctionSpace());
1393     V.typeMatchRight(*this);
1394     m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1395     return make_tuple(boost::python::object(ev),boost::python::object(V));
1396     }
1397    
1398     const boost::python::tuple
1399 jgs 121 Data::mindp() const
1400     {
1401 jgs 148 // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1402     // abort (for unknown reasons) if there are openmp directives with it in the
1403     // surrounding function
1404    
1405     int SampleNo;
1406     int DataPointNo;
1407    
1408     calc_mindp(SampleNo,DataPointNo);
1409    
1410     return make_tuple(SampleNo,DataPointNo);
1411     }
1412    
1413     void
1414     Data::calc_mindp(int& SampleNo,
1415     int& DataPointNo) const
1416     {
1417     int i,j;
1418     int lowi=0,lowj=0;
1419     double min=numeric_limits<double>::max();
1420    
1421 jgs 121 Data temp=minval();
1422    
1423     int numSamples=temp.getNumSamples();
1424     int numDPPSample=temp.getNumDataPointsPerSample();
1425    
1426 jgs 148 double next,local_min;
1427     int local_lowi,local_lowj;
1428 jgs 121
1429 jgs 148 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1430     {
1431     local_min=min;
1432     #pragma omp for private(i,j) schedule(static)
1433     for (i=0; i<numSamples; i++) {
1434     for (j=0; j<numDPPSample; j++) {
1435     next=temp.getDataPoint(i,j)();
1436     if (next<local_min) {
1437     local_min=next;
1438     local_lowi=i;
1439     local_lowj=j;
1440     }
1441 jgs 121 }
1442     }
1443 jgs 148 #pragma omp critical
1444     if (local_min<min) {
1445     min=local_min;
1446     lowi=local_lowi;
1447     lowj=local_lowj;
1448     }
1449 jgs 121 }
1450    
1451 jgs 148 SampleNo = lowi;
1452     DataPointNo = lowj;
1453 jgs 121 }
1454    
1455 jgs 104 void
1456     Data::saveDX(std::string fileName) const
1457     {
1458 jgs 153 boost::python::dict args;
1459     args["data"]=boost::python::object(this);
1460     getDomain().saveDX(fileName,args);
1461 jgs 104 return;
1462     }
1463    
1464 jgs 110 void
1465     Data::saveVTK(std::string fileName) const
1466     {
1467 jgs 153 boost::python::dict args;
1468     args["data"]=boost::python::object(this);
1469     getDomain().saveVTK(fileName,args);
1470 jgs 110 return;
1471     }
1472    
1473 jgs 102 Data&
1474     Data::operator+=(const Data& right)
1475     {
1476 jgs 151 #if defined DOPROF
1477 jgs 123 profData->binary++;
1478 jgs 151 #endif
1479 jgs 94 binaryOp(right,plus<double>());
1480     return (*this);
1481     }
1482    
1483 jgs 102 Data&
1484     Data::operator+=(const boost::python::object& right)
1485 jgs 94 {
1486 jgs 151 #if defined DOPROF
1487 jgs 123 profData->binary++;
1488 jgs 151 #endif
1489 jgs 94 binaryOp(right,plus<double>());
1490     return (*this);
1491     }
1492    
1493 jgs 102 Data&
1494     Data::operator-=(const Data& right)
1495 jgs 94 {
1496 jgs 151 #if defined DOPROF
1497 jgs 123 profData->binary++;
1498 jgs 151 #endif
1499 jgs 94 binaryOp(right,minus<double>());
1500     return (*this);
1501     }
1502    
1503 jgs 102 Data&
1504     Data::operator-=(const boost::python::object& right)
1505 jgs 94 {
1506 jgs 151 #if defined DOPROF
1507 jgs 123 profData->binary++;
1508 jgs 151 #endif
1509 jgs 94 binaryOp(right,minus<double>());
1510     return (*this);
1511     }
1512    
1513 jgs 102 Data&
1514     Data::operator*=(const Data& right)
1515 jgs 94 {
1516 jgs 151 #if defined DOPROF
1517 jgs 123 profData->binary++;
1518 jgs 151 #endif
1519 jgs 94 binaryOp(right,multiplies<double>());
1520     return (*this);
1521     }
1522    
1523 jgs 102 Data&
1524     Data::operator*=(const boost::python::object& right)
1525 jgs 94 {
1526 jgs 151 #if defined DOPROF
1527 jgs 123 profData->binary++;
1528 jgs 151 #endif
1529 jgs 94 binaryOp(right,multiplies<double>());
1530     return (*this);
1531     }
1532    
1533 jgs 102 Data&
1534     Data::operator/=(const Data& right)
1535 jgs 94 {
1536 jgs 151 #if defined DOPROF
1537 jgs 123 profData->binary++;
1538 jgs 151 #endif
1539 jgs 94 binaryOp(right,divides<double>());
1540     return (*this);
1541     }
1542    
1543 jgs 102 Data&
1544     Data::operator/=(const boost::python::object& right)
1545 jgs 94 {
1546 jgs 151 #if defined DOPROF
1547 jgs 123 profData->binary++;
1548 jgs 151 #endif
1549 jgs 94 binaryOp(right,divides<double>());
1550     return (*this);
1551     }
1552    
1553 jgs 102 Data
1554 gross 699 Data::rpowO(const boost::python::object& left) const
1555     {
1556     #if defined DOPROF
1557     profData->binary++;
1558     #endif
1559     Data left_d(left,*this);
1560     return left_d.powD(*this);
1561     }
1562    
1563     Data
1564 jgs 102 Data::powO(const boost::python::object& right) const
1565 jgs 94 {
1566 jgs 151 #if defined DOPROF
1567 jgs 123 profData->binary++;
1568 jgs 151 #endif
1569 jgs 94 Data result;
1570     result.copy(*this);
1571     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1572     return result;
1573     }
1574    
1575 jgs 102 Data
1576     Data::powD(const Data& right) const
1577 jgs 94 {
1578 jgs 151 #if defined DOPROF
1579 jgs 123 profData->binary++;
1580 jgs 151 #endif
1581 jgs 94 Data result;
1582     result.copy(*this);
1583     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1584     return result;
1585     }
1586    
1587 bcumming 751 void
1588     Data::print()
1589     {
1590     int i,j;
1591    
1592     printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
1593     for( i=0; i<getNumSamples(); i++ )
1594     {
1595     printf( "[%6d]", i );
1596     for( j=0; j<getNumDataPointsPerSample(); j++ )
1597     printf( "\t%10.7g", (getSampleData(i))[j] );
1598     printf( "\n" );
1599     }
1600     }
1601    
1602 jgs 94 //
1603 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1604 jgs 102 Data
1605     escript::operator+(const Data& left, const Data& right)
1606 jgs 94 {
1607     Data result;
1608     //
1609     // perform a deep copy
1610     result.copy(left);
1611     result+=right;
1612     return result;
1613     }
1614    
1615     //
1616 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1617 jgs 102 Data
1618     escript::operator-(const Data& left, const Data& right)
1619 jgs 94 {
1620     Data result;
1621     //
1622     // perform a deep copy
1623     result.copy(left);
1624     result-=right;
1625     return result;
1626     }
1627    
1628     //
1629 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1630 jgs 102 Data
1631     escript::operator*(const Data& left, const Data& right)
1632 jgs 94 {
1633     Data result;
1634     //
1635     // perform a deep copy
1636     result.copy(left);
1637     result*=right;
1638     return result;
1639     }
1640    
1641     //
1642 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1643 jgs 102 Data
1644     escript::operator/(const Data& left, const Data& right)
1645 jgs 94 {
1646     Data result;
1647     //
1648     // perform a deep copy
1649     result.copy(left);
1650     result/=right;
1651     return result;
1652     }
1653    
1654     //
1655 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1656 jgs 102 Data
1657     escript::operator+(const Data& left, const boost::python::object& right)
1658 jgs 94 {
1659     //
1660     // Convert to DataArray format if possible
1661     DataArray temp(right);
1662     Data result;
1663     //
1664     // perform a deep copy
1665     result.copy(left);
1666     result+=right;
1667     return result;
1668     }
1669    
1670     //
1671 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1672 jgs 102 Data
1673     escript::operator-(const Data& left, const boost::python::object& right)
1674 jgs 94 {
1675     //
1676     // Convert to DataArray format if possible
1677     DataArray temp(right);
1678     Data result;
1679     //
1680     // perform a deep copy
1681     result.copy(left);
1682     result-=right;
1683     return result;
1684     }
1685    
1686     //
1687 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1688 jgs 102 Data
1689     escript::operator*(const Data& left, const boost::python::object& right)
1690 jgs 94 {
1691     //
1692     // Convert to DataArray format if possible
1693     DataArray temp(right);
1694     Data result;
1695     //
1696     // perform a deep copy
1697     result.copy(left);
1698     result*=right;
1699     return result;
1700     }
1701    
1702     //
1703 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1704 jgs 102 Data
1705     escript::operator/(const Data& left, const boost::python::object& right)
1706 jgs 94 {
1707     //
1708     // Convert to DataArray format if possible
1709     DataArray temp(right);
1710     Data result;
1711     //
1712     // perform a deep copy
1713     result.copy(left);
1714     result/=right;
1715     return result;
1716     }
1717    
1718     //
1719 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1720 jgs 102 Data
1721     escript::operator+(const boost::python::object& left, const Data& right)
1722 jgs 94 {
1723     //
1724     // Construct the result using the given value and the other parameters
1725     // from right
1726     Data result(left,right);
1727     result+=right;
1728     return result;
1729     }
1730    
1731     //
1732 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1733 jgs 102 Data
1734     escript::operator-(const boost::python::object& left, const Data& right)
1735 jgs 94 {
1736     //
1737     // Construct the result using the given value and the other parameters
1738     // from right
1739     Data result(left,right);
1740     result-=right;
1741     return result;
1742     }
1743    
1744     //
1745 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1746 jgs 102 Data
1747     escript::operator*(const boost::python::object& left, const Data& right)
1748 jgs 94 {
1749     //
1750     // Construct the result using the given value and the other parameters
1751     // from right
1752     Data result(left,right);
1753     result*=right;
1754     return result;
1755     }
1756    
1757     //
1758 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1759 jgs 102 Data
1760     escript::operator/(const boost::python::object& left, const Data& right)
1761 jgs 94 {
1762     //
1763     // Construct the result using the given value and the other parameters
1764     // from right
1765     Data result(left,right);
1766     result/=right;
1767     return result;
1768     }
1769    
1770     //
1771 jgs 102 //bool escript::operator==(const Data& left, const Data& right)
1772     //{
1773     // /*
1774     // NB: this operator does very little at this point, and isn't to
1775     // be relied on. Requires further implementation.
1776     // */
1777     //
1778     // bool ret;
1779     //
1780     // if (left.isEmpty()) {
1781     // if(!right.isEmpty()) {
1782     // ret = false;
1783     // } else {
1784     // ret = true;
1785     // }
1786     // }
1787     //
1788     // if (left.isConstant()) {
1789     // if(!right.isConstant()) {
1790     // ret = false;
1791     // } else {
1792     // ret = true;
1793     // }
1794     // }
1795     //
1796     // if (left.isTagged()) {
1797     // if(!right.isTagged()) {
1798     // ret = false;
1799     // } else {
1800     // ret = true;
1801     // }
1802     // }
1803     //
1804     // if (left.isExpanded()) {
1805     // if(!right.isExpanded()) {
1806     // ret = false;
1807     // } else {
1808     // ret = true;
1809     // }
1810     // }
1811     //
1812     // return ret;
1813     //}
1814    
1815 bcumming 751 /* TODO */
1816     /* global reduction */
1817 jgs 102 Data
1818     Data::getItem(const boost::python::object& key) const
1819 jgs 94 {
1820 jgs 102 const DataArrayView& view=getPointDataView();
1821 jgs 94
1822 jgs 102 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1823 jgs 94
1824 jgs 102 if (slice_region.size()!=view.getRank()) {
1825     throw DataException("Error - slice size does not match Data rank.");
1826 jgs 94 }
1827    
1828 jgs 102 return getSlice(slice_region);
1829 jgs 94 }
1830    
1831 bcumming 751 /* TODO */
1832     /* global reduction */
1833 jgs 94 Data
1834 jgs 102 Data::getSlice(const DataArrayView::RegionType& region) const
1835 jgs 94 {
1836 jgs 151 #if defined DOPROF
1837 jgs 123 profData->slicing++;
1838 jgs 151 #endif
1839 jgs 102 return Data(*this,region);
1840 jgs 94 }
1841    
1842 bcumming 751 /* TODO */
1843     /* global reduction */
1844 jgs 94 void
1845 jgs 102 Data::setItemO(const boost::python::object& key,
1846     const boost::python::object& value)
1847 jgs 94 {
1848 jgs 102 Data tempData(value,getFunctionSpace());
1849     setItemD(key,tempData);
1850     }
1851    
1852 bcumming 751 /* TODO */
1853     /* global reduction */
1854 jgs 102 void
1855     Data::setItemD(const boost::python::object& key,
1856     const Data& value)
1857     {
1858 jgs 94 const DataArrayView& view=getPointDataView();
1859 jgs 123
1860 jgs 94 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1861     if (slice_region.size()!=view.getRank()) {
1862     throw DataException("Error - slice size does not match Data rank.");
1863     }
1864 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
1865     setSlice(Data(value,getFunctionSpace()),slice_region);
1866     } else {
1867     setSlice(value,slice_region);
1868     }
1869 jgs 94 }
1870    
1871 bcumming 751 /* TODO */
1872     /* global reduction */
1873 jgs 94 void
1874 jgs 102 Data::setSlice(const Data& value,
1875     const DataArrayView::RegionType& region)
1876 jgs 94 {
1877 jgs 151 #if defined DOPROF
1878 jgs 123 profData->slicing++;
1879 jgs 151 #endif
1880 jgs 102 Data tempValue(value);
1881     typeMatchLeft(tempValue);
1882     typeMatchRight(tempValue);
1883     m_data->setSlice(tempValue.m_data.get(),region);
1884     }
1885    
1886     void
1887     Data::typeMatchLeft(Data& right) const
1888     {
1889     if (isExpanded()){
1890     right.expand();
1891     } else if (isTagged()) {
1892     if (right.isConstant()) {
1893     right.tag();
1894     }
1895     }
1896     }
1897    
1898     void
1899     Data::typeMatchRight(const Data& right)
1900     {
1901 jgs 94 if (isTagged()) {
1902     if (right.isExpanded()) {
1903     expand();
1904     }
1905     } else if (isConstant()) {
1906     if (right.isExpanded()) {
1907     expand();
1908     } else if (right.isTagged()) {
1909     tag();
1910     }
1911     }
1912     }
1913    
1914 bcumming 751 /* TODO */
1915     /* global reduction */
1916 jgs 94 void
1917     Data::setTaggedValue(int tagKey,
1918     const boost::python::object& value)
1919     {
1920     //
1921     // Ensure underlying data object is of type DataTagged
1922     tag();
1923    
1924     if (!isTagged()) {
1925     throw DataException("Error - DataTagged conversion failed!!");
1926     }
1927    
1928     //
1929     // Construct DataArray from boost::python::object input value
1930     DataArray valueDataArray(value);
1931    
1932     //
1933     // Call DataAbstract::setTaggedValue
1934     m_data->setTaggedValue(tagKey,valueDataArray.getView());
1935     }
1936    
1937 bcumming 751 /* TODO */
1938     /* global reduction */
1939 jgs 110 void
1940 jgs 121 Data::setTaggedValueFromCPP(int tagKey,
1941     const DataArrayView& value)
1942     {
1943     //
1944     // Ensure underlying data object is of type DataTagged
1945     tag();
1946    
1947     if (!isTagged()) {
1948     throw DataException("Error - DataTagged conversion failed!!");
1949     }
1950    
1951     //
1952     // Call DataAbstract::setTaggedValue
1953     m_data->setTaggedValue(tagKey,value);
1954     }
1955    
1956 bcumming 751 /* TODO */
1957     /* global reduction */
1958 jgs 149 int
1959     Data::getTagNumber(int dpno)
1960     {
1961     return m_data->getTagNumber(dpno);
1962     }
1963    
1964 bcumming 751 /* TODO */
1965     /* global reduction */
1966 jgs 121 void
1967 jgs 110 Data::setRefValue(int ref,
1968     const boost::python::numeric::array& value)
1969     {
1970     //
1971     // Construct DataArray from boost::python::object input value
1972     DataArray valueDataArray(value);
1973    
1974     //
1975     // Call DataAbstract::setRefValue
1976     m_data->setRefValue(ref,valueDataArray);
1977     }
1978    
1979 bcumming 751 /* TODO */
1980     /* global reduction */
1981 jgs 110 void
1982     Data::getRefValue(int ref,
1983     boost::python::numeric::array& value)
1984     {
1985     //
1986     // Construct DataArray for boost::python::object return value
1987     DataArray valueDataArray(value);
1988    
1989     //
1990     // Load DataArray with values from data-points specified by ref
1991     m_data->getRefValue(ref,valueDataArray);
1992    
1993     //
1994     // Load values from valueDataArray into return numarray
1995    
1996     // extract the shape of the numarray
1997     int rank = value.getrank();
1998     DataArrayView::ShapeType shape;
1999     for (int i=0; i < rank; i++) {
2000     shape.push_back(extract<int>(value.getshape()[i]));
2001     }
2002    
2003     // and load the numarray with the data from the DataArray
2004     DataArrayView valueView = valueDataArray.getView();
2005    
2006     if (rank==0) {
2007 jgs 126 boost::python::numeric::array temp_numArray(valueView());
2008     value = temp_numArray;
2009 jgs 110 }
2010     if (rank==1) {
2011     for (int i=0; i < shape[0]; i++) {
2012     value[i] = valueView(i);
2013     }
2014     }
2015     if (rank==2) {
2016 jgs 126 for (int i=0; i < shape[0]; i++) {
2017     for (int j=0; j < shape[1]; j++) {
2018     value[i][j] = valueView(i,j);
2019     }
2020     }
2021 jgs 110 }
2022     if (rank==3) {
2023 jgs 126 for (int i=0; i < shape[0]; i++) {
2024     for (int j=0; j < shape[1]; j++) {
2025     for (int k=0; k < shape[2]; k++) {
2026     value[i][j][k] = valueView(i,j,k);
2027     }
2028     }
2029     }
2030 jgs 110 }
2031     if (rank==4) {
2032 jgs 126 for (int i=0; i < shape[0]; i++) {
2033     for (int j=0; j < shape[1]; j++) {
2034     for (int k=0; k < shape[2]; k++) {
2035     for (int l=0; l < shape[3]; l++) {
2036     value[i][j][k][l] = valueView(i,j,k,l);
2037     }
2038     }
2039     }
2040     }
2041 jgs 110 }
2042    
2043     }
2044    
2045 jgs 94 void
2046 jgs 119 Data::archiveData(const std::string fileName)
2047     {
2048     cout << "Archiving Data object to: " << fileName << endl;
2049    
2050     //
2051     // Determine type of this Data object
2052     int dataType = -1;
2053    
2054     if (isEmpty()) {
2055     dataType = 0;
2056     cout << "\tdataType: DataEmpty" << endl;
2057     }
2058     if (isConstant()) {
2059     dataType = 1;
2060     cout << "\tdataType: DataConstant" << endl;
2061     }
2062     if (isTagged()) {
2063     dataType = 2;
2064     cout << "\tdataType: DataTagged" << endl;
2065     }
2066     if (isExpanded()) {
2067     dataType = 3;
2068     cout << "\tdataType: DataExpanded" << endl;
2069     }
2070 jgs 123
2071 jgs 119 if (dataType == -1) {
2072     throw DataException("archiveData Error: undefined dataType");
2073     }
2074    
2075     //
2076     // Collect data items common to all Data types
2077     int noSamples = getNumSamples();
2078     int noDPPSample = getNumDataPointsPerSample();
2079     int functionSpaceType = getFunctionSpace().getTypeCode();
2080     int dataPointRank = getDataPointRank();
2081     int dataPointSize = getDataPointSize();
2082     int dataLength = getLength();
2083     DataArrayView::ShapeType dataPointShape = getDataPointShape();
2084 woo409 757 vector<int> referenceNumbers(noSamples);
2085 jgs 119 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2086     referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
2087     }
2088 woo409 757 vector<int> tagNumbers(noSamples);
2089 jgs 119 if (isTagged()) {
2090     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2091     tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
2092     }
2093     }
2094    
2095     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2096     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2097     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2098    
2099     //
2100     // Flatten Shape to an array of integers suitable for writing to file
2101     int flatShape[4] = {0,0,0,0};
2102     cout << "\tshape: < ";
2103     for (int dim=0; dim<dataPointRank; dim++) {
2104     flatShape[dim] = dataPointShape[dim];
2105     cout << dataPointShape[dim] << " ";
2106     }
2107     cout << ">" << endl;
2108    
2109     //
2110 jgs 123 // Open archive file
2111 jgs 119 ofstream archiveFile;
2112     archiveFile.open(fileName.data(), ios::out);
2113    
2114     if (!archiveFile.good()) {
2115     throw DataException("archiveData Error: problem opening archive file");
2116     }
2117    
2118 jgs 123 //
2119     // Write common data items to archive file
2120 jgs 119 archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
2121     archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
2122     archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2123     archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2124     archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2125     archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2126     archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
2127     for (int dim = 0; dim < 4; dim++) {
2128     archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2129     }
2130     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2131     archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2132     }
2133     if (isTagged()) {
2134     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2135     archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2136     }
2137     }
2138    
2139     if (!archiveFile.good()) {
2140     throw DataException("archiveData Error: problem writing to archive file");
2141     }
2142    
2143     //
2144 jgs 123 // Archive underlying data values for each Data type
2145     int noValues;
2146 jgs 119 switch (dataType) {
2147     case 0:
2148     // DataEmpty
2149 jgs 123 noValues = 0;
2150     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2151     cout << "\tnoValues: " << noValues << endl;
2152 jgs 119 break;
2153     case 1:
2154     // DataConstant
2155 jgs 123 noValues = m_data->getLength();
2156     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2157     cout << "\tnoValues: " << noValues << endl;
2158     if (m_data->archiveData(archiveFile,noValues)) {
2159     throw DataException("archiveData Error: problem writing data to archive file");
2160     }
2161 jgs 119 break;
2162     case 2:
2163     // DataTagged
2164 jgs 123 noValues = m_data->getLength();
2165     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2166     cout << "\tnoValues: " << noValues << endl;
2167     if (m_data->archiveData(archiveFile,noValues)) {
2168     throw DataException("archiveData Error: problem writing data to archive file");
2169     }
2170 jgs 119 break;
2171     case 3:
2172     // DataExpanded
2173 jgs 123 noValues = m_data->getLength();
2174     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2175     cout << "\tnoValues: " << noValues << endl;
2176     if (m_data->archiveData(archiveFile,noValues)) {
2177     throw DataException("archiveData Error: problem writing data to archive file");
2178     }
2179 jgs 119 break;
2180     }
2181    
2182 jgs 123 if (!archiveFile.good()) {
2183     throw DataException("archiveData Error: problem writing data to archive file");
2184     }
2185    
2186     //
2187     // Close archive file
2188     archiveFile.close();
2189    
2190     if (!archiveFile.good()) {
2191     throw DataException("archiveData Error: problem closing archive file");
2192     }
2193    
2194 jgs 119 }
2195    
2196     void
2197     Data::extractData(const std::string fileName,
2198     const FunctionSpace& fspace)
2199     {
2200     //
2201     // Can only extract Data to an object which is initially DataEmpty
2202     if (!isEmpty()) {
2203     throw DataException("extractData Error: can only extract to DataEmpty object");
2204     }
2205    
2206     cout << "Extracting Data object from: " << fileName << endl;
2207    
2208     int dataType;
2209     int noSamples;
2210     int noDPPSample;
2211     int functionSpaceType;
2212     int dataPointRank;
2213     int dataPointSize;
2214     int dataLength;
2215     DataArrayView::ShapeType dataPointShape;
2216     int flatShape[4];
2217    
2218     //
2219 jgs 123 // Open the archive file
2220 jgs 119 ifstream archiveFile;
2221     archiveFile.open(fileName.data(), ios::in);
2222    
2223     if (!archiveFile.good()) {
2224     throw DataException("extractData Error: problem opening archive file");
2225     }
2226    
2227 jgs 123 //
2228     // Read common data items from archive file
2229 jgs 119 archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
2230     archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2231     archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2232     archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2233     archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2234     archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2235     archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
2236     for (int dim = 0; dim < 4; dim++) {
2237     archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2238     if (flatShape[dim]>0) {
2239     dataPointShape.push_back(flatShape[dim]);
2240     }
2241     }
2242 woo409 757 vector<int> referenceNumbers(noSamples);
2243 jgs 119 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2244     archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2245     }
2246 woo409 757 vector<int> tagNumbers(noSamples);
2247 jgs 119 if (dataType==2) {
2248     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2249     archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2250     }
2251     }
2252    
2253     if (!archiveFile.good()) {
2254     throw DataException("extractData Error: problem reading from archive file");
2255     }
2256    
2257 jgs 123 //
2258     // Verify the values just read from the archive file
2259 jgs 119 switch (dataType) {
2260     case 0:
2261     cout << "\tdataType: DataEmpty" << endl;
2262     break;
2263     case 1:
2264     cout << "\tdataType: DataConstant" << endl;
2265     break;
2266     case 2:
2267     cout << "\tdataType: DataTagged" << endl;
2268     break;
2269     case 3:
2270     cout << "\tdataType: DataExpanded" << endl;
2271     break;
2272     default:
2273     throw DataException("extractData Error: undefined dataType read from archive file");
2274     break;
2275     }
2276    
2277     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2278     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2279     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2280     cout << "\tshape: < ";
2281     for (int dim = 0; dim < dataPointRank; dim++) {
2282     cout << dataPointShape[dim] << " ";
2283     }
2284     cout << ">" << endl;
2285    
2286     //
2287     // Verify that supplied FunctionSpace object is compatible with this Data object.
2288     if ( (fspace.getTypeCode()!=functionSpaceType) ||
2289     (fspace.getNumSamples()!=noSamples) ||
2290     (fspace.getNumDPPSample()!=noDPPSample)
2291     ) {
2292     throw DataException("extractData Error: incompatible FunctionSpace");
2293     }
2294     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2295     if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2296     throw DataException("extractData Error: incompatible FunctionSpace");
2297     }
2298     }
2299     if (dataType==2) {
2300     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2301     if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2302     throw DataException("extractData Error: incompatible FunctionSpace");
2303     }
2304     }
2305     }
2306    
2307     //
2308     // Construct a DataVector to hold underlying data values
2309     DataVector dataVec(dataLength);
2310    
2311     //
2312     // Load this DataVector with the appropriate values
2313 jgs 123 int noValues;
2314     archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2315     cout << "\tnoValues: " << noValues << endl;
2316 jgs 119 switch (dataType) {
2317     case 0:
2318     // DataEmpty
2319 jgs 123 if (noValues != 0) {
2320     throw DataException("extractData Error: problem reading data from archive file");
2321     }
2322 jgs 119 break;
2323     case 1:
2324     // DataConstant
2325 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2326     throw DataException("extractData Error: problem reading data from archive file");
2327     }
2328 jgs 119 break;
2329     case 2:
2330     // DataTagged
2331 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2332     throw DataException("extractData Error: problem reading data from archive file");
2333     }
2334 jgs 119 break;
2335     case 3:
2336     // DataExpanded
2337 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2338     throw DataException("extractData Error: problem reading data from archive file");
2339     }
2340 jgs 119 break;
2341     }
2342    
2343 jgs 123 if (!archiveFile.good()) {
2344     throw DataException("extractData Error: problem reading from archive file");
2345     }
2346    
2347 jgs 119 //
2348 jgs 123 // Close archive file
2349     archiveFile.close();
2350    
2351     if (!archiveFile.good()) {
2352     throw DataException("extractData Error: problem closing archive file");
2353     }
2354    
2355     //
2356 jgs 119 // Construct an appropriate Data object
2357     DataAbstract* tempData;
2358     switch (dataType) {
2359     case 0:
2360     // DataEmpty
2361     tempData=new DataEmpty();
2362     break;
2363     case 1:
2364     // DataConstant
2365     tempData=new DataConstant(fspace,dataPointShape,dataVec);
2366     break;
2367     case 2:
2368     // DataTagged
2369     tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2370     break;
2371     case 3:
2372     // DataExpanded
2373     tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2374     break;
2375     }
2376     shared_ptr<DataAbstract> temp_data(tempData);
2377     m_data=temp_data;
2378     }
2379    
2380 jgs 94 ostream& escript::operator<<(ostream& o, const Data& data)
2381     {
2382     o << data.toString();
2383     return o;
2384     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26