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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 751 - (hide annotations)
Mon Jun 26 01:46:34 2006 UTC (13 years, 3 months ago) by bcumming
File size: 54964 byte(s)
Changes relating to the MPI version of escript
The standard OpenMP version of escript is unchanged

- updated data types (Finley_Mesh, Finley_NodeFile, etc) to store meshes
  over multiple MPI processes.
- added CommBuffer code in Paso for communication of Data associated
  with distributed meshes
- updates in Finley and Escript to support distributed data and operations
  on distributed data (such as interpolation).
- construction of RHS in MPI, so that simple explicit schemes (such as
  /docs/examples/wave.py without IO and the Locator) can run in MPI.
- updated mesh generation for first order line, rectangle and brick
  meshes and second order line meshes in MPI.        
- small changes to trunk/SConstruct and trunk/scons/ess_options.py to
  build the MPI version, these changes are turned off by default.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26