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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 783 - (hide annotations)
Tue Jul 18 01:32:50 2006 UTC (12 years, 11 months ago) by gross
File size: 62557 byte(s)
coordinates, element size and normals returned by corresponding
FunctionSpace mesthods are now protected against updates. So 
+=, -=, *=, /=, setTaggedValue, fillFromNumArray will through an
excpetion.

The FunctionSpace class does nut buffer the oordinates, element size and
normals yet.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26