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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 790 - (hide annotations)
Wed Jul 26 23:12:34 2006 UTC (13 years, 1 month ago) by bcumming
File size: 64492 byte(s)
changes to escript/py_src/pdetools.py and /escript/src/Data.h/.cpp to
make the Locator work in MPI. escript::Data::mindp now returns a 3 tuple,
with the MPI rank of the process on which the minimum value occurs
included. escript::Data::convertToNumArrayFromDPNo also takes the ProcNo
to perform the MPI reduction.

This had to be implemented in both the MPI and non-MPI versions to allow
the necesary changes to the Python code in pdetools.py. In the non-MPI
version ProcNo is set to 0. This works for the explicit scripts tested
thus far, however if it causes problems in your scripts contact Ben or
Lutz, or revert the three files (pdetools.py, Data.h and Data.cpp) to
the previous version.  


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 Data::convertToNumArrayFromDPNo(int procNo,
786     int sampleNo,
787     int dataPointNo)
788    
789 jgs 121 {
790 bcumming 790 size_t length=0;
791     int i, j, k, l, pos;
792    
793 jgs 121 //
794     // Check a valid sample number has been supplied
795     if (sampleNo >= getNumSamples()) {
796     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
797     }
798    
799     //
800     // Check a valid data point number has been supplied
801     if (dataPointNo >= getNumDataPointsPerSample()) {
802     throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");
803     }
804    
805     //
806     // determine the rank and shape of each data point
807     int dataPointRank = getDataPointRank();
808     DataArrayView::ShapeType dataPointShape = getDataPointShape();
809    
810     //
811     // create the numeric array to be returned
812     boost::python::numeric::array numArray(0.0);
813    
814     //
815     // the shape of the returned numeric array will be the same
816     // as that of the data point
817     int arrayRank = dataPointRank;
818     DataArrayView::ShapeType arrayShape = dataPointShape;
819    
820     //
821     // resize the numeric array to the shape just calculated
822     if (arrayRank==0) {
823     numArray.resize(1);
824     }
825     if (arrayRank==1) {
826     numArray.resize(arrayShape[0]);
827     }
828     if (arrayRank==2) {
829     numArray.resize(arrayShape[0],arrayShape[1]);
830     }
831     if (arrayRank==3) {
832     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
833     }
834     if (arrayRank==4) {
835     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
836     }
837    
838 bcumming 790 // added for the MPI communication
839     length=1;
840     for( i=0; i<arrayRank; i++ )
841     length *= arrayShape[i];
842     double *tmpData = new double[length];
843    
844 jgs 121 //
845     // load the values for the data point into the numeric array.
846 bcumming 790
847     // updated for the MPI case
848     if( get_MPIRank()==procNo ){
849     // create a view of the data if it is stored locally
850     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
851    
852     // pack the data from the view into tmpData for MPI communication
853     pos=0;
854     switch( dataPointRank ){
855     case 0 :
856     tmpData[0] = dataPointView();
857     break;
858     case 1 :
859     for( i=0; i<dataPointShape[0]; i++ )
860     tmpData[i]=dataPointView(i);
861     break;
862     case 2 :
863     for( i=0; i<dataPointShape[0]; i++ )
864     for( j=0; j<dataPointShape[1]; j++, pos++ )
865     tmpData[pos]=dataPointView(i,j);
866     break;
867     case 3 :
868     for( i=0; i<dataPointShape[0]; i++ )
869     for( j=0; j<dataPointShape[1]; j++ )
870     for( k=0; k<dataPointShape[2]; k++, pos++ )
871     tmpData[pos]=dataPointView(i,j,k);
872     break;
873     case 4 :
874     for( i=0; i<dataPointShape[0]; i++ )
875     for( j=0; j<dataPointShape[1]; j++ )
876     for( k=0; k<dataPointShape[2]; k++ )
877     for( l=0; l<dataPointShape[3]; l++, pos++ )
878     tmpData[pos]=dataPointView(i,j,k,l);
879     break;
880     }
881     }
882     #ifdef PASO_MPI
883     // broadcast the data to all other processes
884     MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
885     #endif
886    
887     // unpack the data
888     switch( dataPointRank ){
889     case 0 :
890     numArray[i]=tmpData[0];
891     break;
892     case 1 :
893     for( i=0; i<dataPointShape[0]; i++ )
894     numArray[i]=tmpData[i];
895     break;
896     case 2 :
897     for( i=0; i<dataPointShape[0]; i++ )
898     for( j=0; j<dataPointShape[1]; j++ )
899     tmpData[i+j*dataPointShape[0]];
900     break;
901     case 3 :
902     for( i=0; i<dataPointShape[0]; i++ )
903     for( j=0; j<dataPointShape[1]; j++ )
904     for( k=0; k<dataPointShape[2]; k++ )
905     tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
906     break;
907     case 4 :
908     for( i=0; i<dataPointShape[0]; i++ )
909     for( j=0; j<dataPointShape[1]; j++ )
910     for( k=0; k<dataPointShape[2]; k++ )
911     for( l=0; l<dataPointShape[3]; l++ )
912     tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
913     break;
914     }
915    
916     delete [] tmpData;
917     /*
918 jgs 121 if (dataPointRank==0) {
919     numArray[0]=dataPointView();
920     }
921     if (dataPointRank==1) {
922     for (int i=0; i<dataPointShape[0]; i++) {
923     numArray[i]=dataPointView(i);
924     }
925     }
926     if (dataPointRank==2) {
927     for (int i=0; i<dataPointShape[0]; i++) {
928     for (int j=0; j<dataPointShape[1]; j++) {
929     numArray[i][j] = dataPointView(i,j);
930     }
931     }
932     }
933     if (dataPointRank==3) {
934     for (int i=0; i<dataPointShape[0]; i++) {
935     for (int j=0; j<dataPointShape[1]; j++) {
936     for (int k=0; k<dataPointShape[2]; k++) {
937     numArray[i][j][k]=dataPointView(i,j,k);
938     }
939     }
940     }
941     }
942     if (dataPointRank==4) {
943     for (int i=0; i<dataPointShape[0]; i++) {
944     for (int j=0; j<dataPointShape[1]; j++) {
945     for (int k=0; k<dataPointShape[2]; k++) {
946     for (int l=0; l<dataPointShape[3]; l++) {
947     numArray[i][j][k][l]=dataPointView(i,j,k,l);
948     }
949     }
950     }
951     }
952     }
953 bcumming 790 */
954 jgs 121
955     //
956     // return the loaded array
957     return numArray;
958     }
959    
960     boost::python::numeric::array
961 jgs 94 Data::integrate() const
962     {
963     int index;
964     int rank = getDataPointRank();
965     DataArrayView::ShapeType shape = getDataPointShape();
966    
967 jgs 151 #if defined DOPROF
968 jgs 123 profData->integrate++;
969 jgs 151 #endif
970 jgs 123
971 jgs 94 //
972     // calculate the integral values
973     vector<double> integrals(getDataPointSize());
974     AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);
975    
976     //
977     // create the numeric array to be returned
978     // and load the array with the integral values
979     boost::python::numeric::array bp_array(1.0);
980     if (rank==0) {
981 jgs 108 bp_array.resize(1);
982 jgs 94 index = 0;
983     bp_array[0] = integrals[index];
984     }
985     if (rank==1) {
986     bp_array.resize(shape[0]);
987     for (int i=0; i<shape[0]; i++) {
988     index = i;
989     bp_array[i] = integrals[index];
990     }
991     }
992     if (rank==2) {
993 gross 436 bp_array.resize(shape[0],shape[1]);
994     for (int i=0; i<shape[0]; i++) {
995     for (int j=0; j<shape[1]; j++) {
996     index = i + shape[0] * j;
997     bp_array[make_tuple(i,j)] = integrals[index];
998     }
999     }
1000 jgs 94 }
1001     if (rank==3) {
1002     bp_array.resize(shape[0],shape[1],shape[2]);
1003     for (int i=0; i<shape[0]; i++) {
1004     for (int j=0; j<shape[1]; j++) {
1005     for (int k=0; k<shape[2]; k++) {
1006     index = i + shape[0] * ( j + shape[1] * k );
1007 gross 436 bp_array[make_tuple(i,j,k)] = integrals[index];
1008 jgs 94 }
1009     }
1010     }
1011     }
1012     if (rank==4) {
1013     bp_array.resize(shape[0],shape[1],shape[2],shape[3]);
1014     for (int i=0; i<shape[0]; i++) {
1015     for (int j=0; j<shape[1]; j++) {
1016     for (int k=0; k<shape[2]; k++) {
1017     for (int l=0; l<shape[3]; l++) {
1018     index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
1019 gross 436 bp_array[make_tuple(i,j,k,l)] = integrals[index];
1020 jgs 94 }
1021     }
1022     }
1023     }
1024     }
1025    
1026     //
1027     // return the loaded array
1028     return bp_array;
1029     }
1030    
1031     Data
1032     Data::sin() const
1033     {
1034 jgs 151 #if defined DOPROF
1035 jgs 123 profData->unary++;
1036 jgs 151 #endif
1037 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
1038     }
1039    
1040     Data
1041     Data::cos() const
1042     {
1043 jgs 151 #if defined DOPROF
1044 jgs 123 profData->unary++;
1045 jgs 151 #endif
1046 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
1047     }
1048    
1049     Data
1050     Data::tan() const
1051     {
1052 jgs 151 #if defined DOPROF
1053 jgs 123 profData->unary++;
1054 jgs 151 #endif
1055 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
1056     }
1057    
1058     Data
1059 jgs 150 Data::asin() const
1060     {
1061 jgs 151 #if defined DOPROF
1062 jgs 150 profData->unary++;
1063 jgs 151 #endif
1064 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
1065     }
1066    
1067     Data
1068     Data::acos() const
1069     {
1070 jgs 151 #if defined DOPROF
1071 jgs 150 profData->unary++;
1072 jgs 151 #endif
1073 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
1074     }
1075    
1076     Data
1077     Data::atan() const
1078     {
1079 jgs 151 #if defined DOPROF
1080 jgs 150 profData->unary++;
1081 jgs 151 #endif
1082 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
1083     }
1084    
1085     Data
1086     Data::sinh() const
1087     {
1088 jgs 151 #if defined DOPROF
1089 jgs 150 profData->unary++;
1090 jgs 151 #endif
1091 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
1092     }
1093    
1094     Data
1095     Data::cosh() const
1096     {
1097 jgs 151 #if defined DOPROF
1098 jgs 150 profData->unary++;
1099 jgs 151 #endif
1100 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
1101     }
1102    
1103     Data
1104     Data::tanh() const
1105     {
1106 jgs 151 #if defined DOPROF
1107 jgs 150 profData->unary++;
1108 jgs 151 #endif
1109 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
1110     }
1111    
1112     Data
1113     Data::asinh() const
1114     {
1115 jgs 151 #if defined DOPROF
1116 jgs 150 profData->unary++;
1117 jgs 151 #endif
1118 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
1119     }
1120    
1121     Data
1122     Data::acosh() const
1123     {
1124 jgs 151 #if defined DOPROF
1125 jgs 150 profData->unary++;
1126 jgs 151 #endif
1127 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
1128     }
1129    
1130     Data
1131     Data::atanh() const
1132     {
1133 jgs 151 #if defined DOPROF
1134 jgs 150 profData->unary++;
1135 jgs 151 #endif
1136 jgs 150 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
1137     }
1138    
1139     Data
1140 gross 286 Data::log10() const
1141 jgs 94 {
1142 jgs 151 #if defined DOPROF
1143 jgs 123 profData->unary++;
1144 jgs 151 #endif
1145 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
1146     }
1147    
1148     Data
1149 gross 286 Data::log() const
1150 jgs 94 {
1151 jgs 151 #if defined DOPROF
1152 jgs 123 profData->unary++;
1153 jgs 151 #endif
1154 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
1155     }
1156    
1157 jgs 106 Data
1158     Data::sign() const
1159 jgs 94 {
1160 jgs 151 #if defined DOPROF
1161 jgs 123 profData->unary++;
1162 jgs 151 #endif
1163 jgs 106 return escript::unaryOp(*this,escript::fsign);
1164 jgs 94 }
1165    
1166 jgs 106 Data
1167     Data::abs() const
1168 jgs 94 {
1169 jgs 151 #if defined DOPROF
1170 jgs 123 profData->unary++;
1171 jgs 151 #endif
1172 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
1173 jgs 94 }
1174    
1175 jgs 106 Data
1176     Data::neg() const
1177 jgs 94 {
1178 jgs 151 #if defined DOPROF
1179 jgs 123 profData->unary++;
1180 jgs 151 #endif
1181 jgs 106 return escript::unaryOp(*this,negate<double>());
1182 jgs 94 }
1183    
1184 jgs 102 Data
1185 jgs 106 Data::pos() const
1186 jgs 94 {
1187 jgs 151 #if defined DOPROF
1188 jgs 123 profData->unary++;
1189 jgs 151 #endif
1190 jgs 148 Data result;
1191     // perform a deep copy
1192     result.copy(*this);
1193     return result;
1194 jgs 102 }
1195    
1196     Data
1197 jgs 106 Data::exp() const
1198 jgs 102 {
1199 jgs 151 #if defined DOPROF
1200 jgs 123 profData->unary++;
1201 jgs 151 #endif
1202 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1203 jgs 102 }
1204    
1205     Data
1206 jgs 106 Data::sqrt() const
1207 jgs 102 {
1208 jgs 151 #if defined DOPROF
1209 jgs 123 profData->unary++;
1210 jgs 151 #endif
1211 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1212 jgs 102 }
1213    
1214 jgs 106 double
1215     Data::Lsup() const
1216 jgs 102 {
1217 bcumming 751 double localValue, globalValue;
1218 jgs 151 #if defined DOPROF
1219 jgs 123 profData->reduction1++;
1220 jgs 151 #endif
1221 jgs 106 //
1222     // set the initial absolute maximum value to zero
1223 bcumming 751
1224 jgs 147 AbsMax abs_max_func;
1225 bcumming 751 localValue = algorithm(abs_max_func,0);
1226     #ifdef PASO_MPI
1227     MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1228     return globalValue;
1229     #else
1230     return localValue;
1231     #endif
1232 jgs 102 }
1233    
1234 jgs 106 double
1235 jgs 117 Data::Linf() const
1236     {
1237 bcumming 751 double localValue, globalValue;
1238 jgs 151 #if defined DOPROF
1239 jgs 123 profData->reduction1++;
1240 jgs 151 #endif
1241 jgs 117 //
1242     // set the initial absolute minimum value to max double
1243 jgs 147 AbsMin abs_min_func;
1244 bcumming 751 localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1245    
1246     #ifdef PASO_MPI
1247     MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1248     return globalValue;
1249     #else
1250     return localValue;
1251     #endif
1252 jgs 117 }
1253    
1254     double
1255 jgs 106 Data::sup() const
1256 jgs 102 {
1257 bcumming 751 double localValue, globalValue;
1258 jgs 151 #if defined DOPROF
1259 jgs 123 profData->reduction1++;
1260 jgs 151 #endif
1261 jgs 106 //
1262     // set the initial maximum value to min possible double
1263 jgs 147 FMax fmax_func;
1264 bcumming 751 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1265     #ifdef PASO_MPI
1266     MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1267     return globalValue;
1268     #else
1269     return localValue;
1270     #endif
1271 jgs 102 }
1272    
1273 jgs 106 double
1274     Data::inf() const
1275 jgs 102 {
1276 bcumming 751 double localValue, globalValue;
1277 jgs 151 #if defined DOPROF
1278 jgs 123 profData->reduction1++;
1279 jgs 151 #endif
1280 jgs 106 //
1281     // set the initial minimum value to max possible double
1282 jgs 147 FMin fmin_func;
1283 bcumming 751 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1284     #ifdef PASO_MPI
1285     MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1286     return globalValue;
1287     #else
1288     return localValue;
1289     #endif
1290 jgs 102 }
1291    
1292 bcumming 751 /* TODO */
1293     /* global reduction */
1294 jgs 102 Data
1295 jgs 106 Data::maxval() const
1296 jgs 102 {
1297 jgs 151 #if defined DOPROF
1298 jgs 123 profData->reduction2++;
1299 jgs 151 #endif
1300 jgs 113 //
1301     // set the initial maximum value to min possible double
1302 jgs 147 FMax fmax_func;
1303     return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1304 jgs 102 }
1305    
1306     Data
1307 jgs 106 Data::minval() const
1308 jgs 102 {
1309 jgs 151 #if defined DOPROF
1310 jgs 123 profData->reduction2++;
1311 jgs 151 #endif
1312 jgs 113 //
1313     // set the initial minimum value to max possible double
1314 jgs 147 FMin fmin_func;
1315     return dp_algorithm(fmin_func,numeric_limits<double>::max());
1316 jgs 102 }
1317    
1318 jgs 123 Data
1319     Data::trace() const
1320     {
1321 jgs 151 #if defined DOPROF
1322 jgs 123 profData->reduction2++;
1323 jgs 151 #endif
1324 jgs 147 Trace trace_func;
1325     return dp_algorithm(trace_func,0);
1326 jgs 123 }
1327    
1328     Data
1329 ksteube 775 Data::symmetric() const
1330 jgs 123 {
1331 ksteube 775 #if defined DOPROF
1332     profData->unary++;
1333     #endif
1334     // check input
1335     DataArrayView::ShapeType s=getDataPointShape();
1336     if (getDataPointRank()==2) {
1337     if(s[0] != s[1])
1338     throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1339     }
1340     else if (getDataPointRank()==4) {
1341     if(!(s[0] == s[2] && s[1] == s[3]))
1342     throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1343     }
1344     else {
1345     throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1346     }
1347     Data ev(0.,getDataPointShape(),getFunctionSpace());
1348     ev.typeMatchRight(*this);
1349     m_data->symmetric(ev.m_data.get());
1350     return ev;
1351     }
1352    
1353     Data
1354     Data::nonsymmetric() const
1355     {
1356     #if defined DOPROF
1357     profData->unary++;
1358     #endif
1359     // check input
1360     DataArrayView::ShapeType s=getDataPointShape();
1361     if (getDataPointRank()==2) {
1362     if(s[0] != s[1])
1363     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1364     DataArrayView::ShapeType ev_shape;
1365     ev_shape.push_back(s[0]);
1366     ev_shape.push_back(s[1]);
1367     Data ev(0.,ev_shape,getFunctionSpace());
1368     ev.typeMatchRight(*this);
1369     m_data->nonsymmetric(ev.m_data.get());
1370     return ev;
1371     }
1372     else if (getDataPointRank()==4) {
1373     if(!(s[0] == s[2] && s[1] == s[3]))
1374     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1375     DataArrayView::ShapeType ev_shape;
1376     ev_shape.push_back(s[0]);
1377     ev_shape.push_back(s[1]);
1378     ev_shape.push_back(s[2]);
1379     ev_shape.push_back(s[3]);
1380     Data ev(0.,ev_shape,getFunctionSpace());
1381     ev.typeMatchRight(*this);
1382     m_data->nonsymmetric(ev.m_data.get());
1383     return ev;
1384     }
1385     else {
1386     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1387     }
1388     }
1389    
1390     Data
1391     Data::matrixtrace(int axis_offset) const
1392     {
1393     #if defined DOPROF
1394     profData->unary++;
1395     #endif
1396     DataArrayView::ShapeType s=getDataPointShape();
1397     if (getDataPointRank()==2) {
1398     DataArrayView::ShapeType ev_shape;
1399     Data ev(0.,ev_shape,getFunctionSpace());
1400     ev.typeMatchRight(*this);
1401     m_data->matrixtrace(ev.m_data.get(), axis_offset);
1402     return ev;
1403     }
1404     if (getDataPointRank()==3) {
1405     DataArrayView::ShapeType ev_shape;
1406     if (axis_offset==0) {
1407     int s2=s[2];
1408     ev_shape.push_back(s2);
1409     }
1410     else if (axis_offset==1) {
1411     int s0=s[0];
1412     ev_shape.push_back(s0);
1413     }
1414     Data ev(0.,ev_shape,getFunctionSpace());
1415     ev.typeMatchRight(*this);
1416     m_data->matrixtrace(ev.m_data.get(), axis_offset);
1417     return ev;
1418     }
1419     if (getDataPointRank()==4) {
1420     DataArrayView::ShapeType ev_shape;
1421     if (axis_offset==0) {
1422     ev_shape.push_back(s[2]);
1423     ev_shape.push_back(s[3]);
1424     }
1425     else if (axis_offset==1) {
1426     ev_shape.push_back(s[0]);
1427     ev_shape.push_back(s[3]);
1428     }
1429     else if (axis_offset==2) {
1430     ev_shape.push_back(s[0]);
1431     ev_shape.push_back(s[1]);
1432     }
1433     Data ev(0.,ev_shape,getFunctionSpace());
1434     ev.typeMatchRight(*this);
1435     m_data->matrixtrace(ev.m_data.get(), axis_offset);
1436     return ev;
1437     }
1438     else {
1439     throw DataException("Error - Data::matrixtrace can only be calculated for rank 2, 3 or 4 object.");
1440     }
1441     }
1442    
1443     Data
1444     Data::transpose(int axis_offset) const
1445     {
1446 jgs 151 #if defined DOPROF
1447 ksteube 775 profData->reduction2++;
1448 jgs 151 #endif
1449 ksteube 775 DataArrayView::ShapeType s=getDataPointShape();
1450     DataArrayView::ShapeType ev_shape;
1451     // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1452     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1453     int rank=getDataPointRank();
1454     if (axis_offset<0 || axis_offset>rank) {
1455     throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1456     }
1457     for (int i=0; i<rank; i++) {
1458     int index = (axis_offset+i)%rank;
1459     ev_shape.push_back(s[index]); // Append to new shape
1460     }
1461     Data ev(0.,ev_shape,getFunctionSpace());
1462     ev.typeMatchRight(*this);
1463     m_data->transpose(ev.m_data.get(), axis_offset);
1464     return ev;
1465 jgs 123 }
1466    
1467 gross 576 Data
1468     Data::eigenvalues() const
1469     {
1470     #if defined DOPROF
1471     profData->unary++;
1472     #endif
1473     // check input
1474     DataArrayView::ShapeType s=getDataPointShape();
1475     if (getDataPointRank()!=2)
1476     throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1477     if(s[0] != s[1])
1478     throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1479     // create return
1480     DataArrayView::ShapeType ev_shape(1,s[0]);
1481     Data ev(0.,ev_shape,getFunctionSpace());
1482     ev.typeMatchRight(*this);
1483     m_data->eigenvalues(ev.m_data.get());
1484     return ev;
1485     }
1486    
1487 jgs 121 const boost::python::tuple
1488 gross 576 Data::eigenvalues_and_eigenvectors(const double tol) const
1489     {
1490     #if defined DOPROF
1491     profData->unary++;
1492     #endif
1493     DataArrayView::ShapeType s=getDataPointShape();
1494     if (getDataPointRank()!=2)
1495     throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1496     if(s[0] != s[1])
1497     throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1498     // create return
1499     DataArrayView::ShapeType ev_shape(1,s[0]);
1500     Data ev(0.,ev_shape,getFunctionSpace());
1501     ev.typeMatchRight(*this);
1502     DataArrayView::ShapeType V_shape(2,s[0]);
1503     Data V(0.,V_shape,getFunctionSpace());
1504     V.typeMatchRight(*this);
1505     m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1506     return make_tuple(boost::python::object(ev),boost::python::object(V));
1507     }
1508    
1509     const boost::python::tuple
1510 jgs 121 Data::mindp() const
1511     {
1512 jgs 148 // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1513     // abort (for unknown reasons) if there are openmp directives with it in the
1514     // surrounding function
1515    
1516     int SampleNo;
1517     int DataPointNo;
1518 bcumming 782 int ProcNo;
1519     calc_mindp(ProcNo,SampleNo,DataPointNo);
1520     return make_tuple(ProcNo,SampleNo,DataPointNo);
1521 jgs 148 }
1522    
1523     void
1524 bcumming 782 Data::calc_mindp( int& ProcNo,
1525     int& SampleNo,
1526     int& DataPointNo) const
1527 jgs 148 {
1528     int i,j;
1529     int lowi=0,lowj=0;
1530     double min=numeric_limits<double>::max();
1531    
1532 jgs 121 Data temp=minval();
1533    
1534     int numSamples=temp.getNumSamples();
1535     int numDPPSample=temp.getNumDataPointsPerSample();
1536    
1537 jgs 148 double next,local_min;
1538     int local_lowi,local_lowj;
1539 jgs 121
1540 jgs 148 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1541     {
1542     local_min=min;
1543     #pragma omp for private(i,j) schedule(static)
1544     for (i=0; i<numSamples; i++) {
1545     for (j=0; j<numDPPSample; j++) {
1546     next=temp.getDataPoint(i,j)();
1547     if (next<local_min) {
1548     local_min=next;
1549     local_lowi=i;
1550     local_lowj=j;
1551     }
1552 jgs 121 }
1553     }
1554 jgs 148 #pragma omp critical
1555     if (local_min<min) {
1556     min=local_min;
1557     lowi=local_lowi;
1558     lowj=local_lowj;
1559     }
1560 jgs 121 }
1561    
1562 bcumming 782 #ifdef PASO_MPI
1563     // determine the processor on which the minimum occurs
1564     next = temp.getDataPoint(lowi,lowj)();
1565     int lowProc = 0;
1566     double *globalMins = new double[get_MPISize()+1];
1567     int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1568    
1569     if( get_MPIRank()==0 ){
1570     next = globalMins[lowProc];
1571     for( i=1; i<get_MPISize(); i++ )
1572     if( next>globalMins[i] ){
1573     lowProc = i;
1574     next = globalMins[i];
1575     }
1576     }
1577     MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1578    
1579     delete [] globalMins;
1580     ProcNo = lowProc;
1581 bcumming 790 #else
1582     ProcNo = 0;
1583 bcumming 782 #endif
1584 jgs 148 SampleNo = lowi;
1585     DataPointNo = lowj;
1586 jgs 121 }
1587    
1588 jgs 104 void
1589     Data::saveDX(std::string fileName) const
1590     {
1591 jgs 153 boost::python::dict args;
1592     args["data"]=boost::python::object(this);
1593     getDomain().saveDX(fileName,args);
1594 jgs 104 return;
1595     }
1596    
1597 jgs 110 void
1598     Data::saveVTK(std::string fileName) const
1599     {
1600 jgs 153 boost::python::dict args;
1601     args["data"]=boost::python::object(this);
1602     getDomain().saveVTK(fileName,args);
1603 jgs 110 return;
1604     }
1605    
1606 jgs 102 Data&
1607     Data::operator+=(const Data& right)
1608     {
1609 gross 783 if (isProtected()) {
1610     throw DataException("Error - attempt to update protected Data object.");
1611     }
1612 jgs 151 #if defined DOPROF
1613 jgs 123 profData->binary++;
1614 jgs 151 #endif
1615 jgs 94 binaryOp(right,plus<double>());
1616     return (*this);
1617     }
1618    
1619 jgs 102 Data&
1620     Data::operator+=(const boost::python::object& right)
1621 jgs 94 {
1622 gross 783 if (isProtected()) {
1623     throw DataException("Error - attempt to update protected Data object.");
1624     }
1625 jgs 151 #if defined DOPROF
1626 jgs 123 profData->binary++;
1627 jgs 151 #endif
1628 jgs 94 binaryOp(right,plus<double>());
1629     return (*this);
1630     }
1631    
1632 jgs 102 Data&
1633     Data::operator-=(const Data& right)
1634 jgs 94 {
1635 gross 783 if (isProtected()) {
1636     throw DataException("Error - attempt to update protected Data object.");
1637     }
1638 jgs 151 #if defined DOPROF
1639 jgs 123 profData->binary++;
1640 jgs 151 #endif
1641 jgs 94 binaryOp(right,minus<double>());
1642     return (*this);
1643     }
1644    
1645 jgs 102 Data&
1646     Data::operator-=(const boost::python::object& right)
1647 jgs 94 {
1648 gross 783 if (isProtected()) {
1649     throw DataException("Error - attempt to update protected Data object.");
1650     }
1651 jgs 151 #if defined DOPROF
1652 jgs 123 profData->binary++;
1653 jgs 151 #endif
1654 jgs 94 binaryOp(right,minus<double>());
1655     return (*this);
1656     }
1657    
1658 jgs 102 Data&
1659     Data::operator*=(const Data& right)
1660 jgs 94 {
1661 gross 783 if (isProtected()) {
1662     throw DataException("Error - attempt to update protected Data object.");
1663     }
1664 jgs 151 #if defined DOPROF
1665 jgs 123 profData->binary++;
1666 jgs 151 #endif
1667 jgs 94 binaryOp(right,multiplies<double>());
1668     return (*this);
1669     }
1670    
1671 jgs 102 Data&
1672     Data::operator*=(const boost::python::object& right)
1673 jgs 94 {
1674 gross 783 if (isProtected()) {
1675     throw DataException("Error - attempt to update protected Data object.");
1676     }
1677 jgs 151 #if defined DOPROF
1678 jgs 123 profData->binary++;
1679 jgs 151 #endif
1680 jgs 94 binaryOp(right,multiplies<double>());
1681     return (*this);
1682     }
1683    
1684 jgs 102 Data&
1685     Data::operator/=(const Data& right)
1686 jgs 94 {
1687 gross 783 if (isProtected()) {
1688     throw DataException("Error - attempt to update protected Data object.");
1689     }
1690 jgs 151 #if defined DOPROF
1691 jgs 123 profData->binary++;
1692 jgs 151 #endif
1693 jgs 94 binaryOp(right,divides<double>());
1694     return (*this);
1695     }
1696    
1697 jgs 102 Data&
1698     Data::operator/=(const boost::python::object& right)
1699 jgs 94 {
1700 gross 783 if (isProtected()) {
1701     throw DataException("Error - attempt to update protected Data object.");
1702     }
1703 jgs 151 #if defined DOPROF
1704 jgs 123 profData->binary++;
1705 jgs 151 #endif
1706 jgs 94 binaryOp(right,divides<double>());
1707     return (*this);
1708     }
1709    
1710 jgs 102 Data
1711 gross 699 Data::rpowO(const boost::python::object& left) const
1712     {
1713 gross 783 if (isProtected()) {
1714     throw DataException("Error - attempt to update protected Data object.");
1715     }
1716 gross 699 #if defined DOPROF
1717     profData->binary++;
1718     #endif
1719     Data left_d(left,*this);
1720     return left_d.powD(*this);
1721     }
1722    
1723     Data
1724 jgs 102 Data::powO(const boost::python::object& right) const
1725 jgs 94 {
1726 jgs 151 #if defined DOPROF
1727 jgs 123 profData->binary++;
1728 jgs 151 #endif
1729 jgs 94 Data result;
1730     result.copy(*this);
1731     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1732     return result;
1733     }
1734    
1735 jgs 102 Data
1736     Data::powD(const Data& right) const
1737 jgs 94 {
1738 jgs 151 #if defined DOPROF
1739 jgs 123 profData->binary++;
1740 jgs 151 #endif
1741 jgs 94 Data result;
1742     result.copy(*this);
1743     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1744     return result;
1745     }
1746    
1747 bcumming 751
1748 jgs 94 //
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 Data& right)
1752 jgs 94 {
1753     Data result;
1754     //
1755     // perform a deep copy
1756     result.copy(left);
1757     result+=right;
1758     return result;
1759     }
1760    
1761     //
1762 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1763 jgs 102 Data
1764     escript::operator-(const Data& left, const Data& right)
1765 jgs 94 {
1766     Data result;
1767     //
1768     // perform a deep copy
1769     result.copy(left);
1770     result-=right;
1771     return result;
1772     }
1773    
1774     //
1775 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1776 jgs 102 Data
1777     escript::operator*(const Data& left, const Data& right)
1778 jgs 94 {
1779     Data result;
1780     //
1781     // perform a deep copy
1782     result.copy(left);
1783     result*=right;
1784     return result;
1785     }
1786    
1787     //
1788 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1789 jgs 102 Data
1790     escript::operator/(const Data& left, const Data& right)
1791 jgs 94 {
1792     Data result;
1793     //
1794     // perform a deep copy
1795     result.copy(left);
1796     result/=right;
1797     return result;
1798     }
1799    
1800     //
1801 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1802 jgs 102 Data
1803     escript::operator+(const Data& left, const boost::python::object& right)
1804 jgs 94 {
1805     //
1806     // Convert to DataArray format if possible
1807     DataArray temp(right);
1808     Data result;
1809     //
1810     // perform a deep copy
1811     result.copy(left);
1812     result+=right;
1813     return result;
1814     }
1815    
1816     //
1817 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1818 jgs 102 Data
1819     escript::operator-(const Data& left, const boost::python::object& right)
1820 jgs 94 {
1821     //
1822     // Convert to DataArray format if possible
1823     DataArray temp(right);
1824     Data result;
1825     //
1826     // perform a deep copy
1827     result.copy(left);
1828     result-=right;
1829     return result;
1830     }
1831    
1832     //
1833 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1834 jgs 102 Data
1835     escript::operator*(const Data& left, const boost::python::object& right)
1836 jgs 94 {
1837     //
1838     // Convert to DataArray format if possible
1839     DataArray temp(right);
1840     Data result;
1841     //
1842     // perform a deep copy
1843     result.copy(left);
1844     result*=right;
1845     return result;
1846     }
1847    
1848     //
1849 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1850 jgs 102 Data
1851     escript::operator/(const Data& left, const boost::python::object& right)
1852 jgs 94 {
1853     //
1854     // Convert to DataArray format if possible
1855     DataArray temp(right);
1856     Data result;
1857     //
1858     // perform a deep copy
1859     result.copy(left);
1860     result/=right;
1861     return result;
1862     }
1863    
1864     //
1865 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1866 jgs 102 Data
1867     escript::operator+(const boost::python::object& left, const Data& right)
1868 jgs 94 {
1869     //
1870     // Construct the result using the given value and the other parameters
1871     // from right
1872     Data result(left,right);
1873     result+=right;
1874     return result;
1875     }
1876    
1877     //
1878 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1879 jgs 102 Data
1880     escript::operator-(const boost::python::object& left, const Data& right)
1881 jgs 94 {
1882     //
1883     // Construct the result using the given value and the other parameters
1884     // from right
1885     Data result(left,right);
1886     result-=right;
1887     return result;
1888     }
1889    
1890     //
1891 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1892 jgs 102 Data
1893     escript::operator*(const boost::python::object& left, const Data& right)
1894 jgs 94 {
1895     //
1896     // Construct the result using the given value and the other parameters
1897     // from right
1898     Data result(left,right);
1899     result*=right;
1900     return result;
1901     }
1902    
1903     //
1904 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1905 jgs 102 Data
1906     escript::operator/(const boost::python::object& left, const Data& right)
1907 jgs 94 {
1908     //
1909     // Construct the result using the given value and the other parameters
1910     // from right
1911     Data result(left,right);
1912     result/=right;
1913     return result;
1914     }
1915    
1916     //
1917 jgs 102 //bool escript::operator==(const Data& left, const Data& right)
1918     //{
1919     // /*
1920     // NB: this operator does very little at this point, and isn't to
1921     // be relied on. Requires further implementation.
1922     // */
1923     //
1924     // bool ret;
1925     //
1926     // if (left.isEmpty()) {
1927     // if(!right.isEmpty()) {
1928     // ret = false;
1929     // } else {
1930     // ret = true;
1931     // }
1932     // }
1933     //
1934     // if (left.isConstant()) {
1935     // if(!right.isConstant()) {
1936     // ret = false;
1937     // } else {
1938     // ret = true;
1939     // }
1940     // }
1941     //
1942     // if (left.isTagged()) {
1943     // if(!right.isTagged()) {
1944     // ret = false;
1945     // } else {
1946     // ret = true;
1947     // }
1948     // }
1949     //
1950     // if (left.isExpanded()) {
1951     // if(!right.isExpanded()) {
1952     // ret = false;
1953     // } else {
1954     // ret = true;
1955     // }
1956     // }
1957     //
1958     // return ret;
1959     //}
1960    
1961 bcumming 751 /* TODO */
1962     /* global reduction */
1963 jgs 102 Data
1964     Data::getItem(const boost::python::object& key) const
1965 jgs 94 {
1966 jgs 102 const DataArrayView& view=getPointDataView();
1967 jgs 94
1968 jgs 102 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1969 jgs 94
1970 jgs 102 if (slice_region.size()!=view.getRank()) {
1971     throw DataException("Error - slice size does not match Data rank.");
1972 jgs 94 }
1973    
1974 jgs 102 return getSlice(slice_region);
1975 jgs 94 }
1976    
1977 bcumming 751 /* TODO */
1978     /* global reduction */
1979 jgs 94 Data
1980 jgs 102 Data::getSlice(const DataArrayView::RegionType& region) const
1981 jgs 94 {
1982 jgs 151 #if defined DOPROF
1983 jgs 123 profData->slicing++;
1984 jgs 151 #endif
1985 jgs 102 return Data(*this,region);
1986 jgs 94 }
1987    
1988 bcumming 751 /* TODO */
1989     /* global reduction */
1990 jgs 94 void
1991 jgs 102 Data::setItemO(const boost::python::object& key,
1992     const boost::python::object& value)
1993 jgs 94 {
1994 jgs 102 Data tempData(value,getFunctionSpace());
1995     setItemD(key,tempData);
1996     }
1997    
1998 bcumming 751 /* TODO */
1999     /* global reduction */
2000 jgs 102 void
2001     Data::setItemD(const boost::python::object& key,
2002     const Data& value)
2003     {
2004 jgs 94 const DataArrayView& view=getPointDataView();
2005 jgs 123
2006 jgs 94 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
2007     if (slice_region.size()!=view.getRank()) {
2008     throw DataException("Error - slice size does not match Data rank.");
2009     }
2010 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
2011     setSlice(Data(value,getFunctionSpace()),slice_region);
2012     } else {
2013     setSlice(value,slice_region);
2014     }
2015 jgs 94 }
2016    
2017 bcumming 751 /* TODO */
2018     /* global reduction */
2019 jgs 94 void
2020 jgs 102 Data::setSlice(const Data& value,
2021     const DataArrayView::RegionType& region)
2022 jgs 94 {
2023 gross 783 if (isProtected()) {
2024     throw DataException("Error - attempt to update protected Data object.");
2025     }
2026 jgs 151 #if defined DOPROF
2027 jgs 123 profData->slicing++;
2028 jgs 151 #endif
2029 jgs 102 Data tempValue(value);
2030     typeMatchLeft(tempValue);
2031     typeMatchRight(tempValue);
2032     m_data->setSlice(tempValue.m_data.get(),region);
2033     }
2034    
2035     void
2036     Data::typeMatchLeft(Data& right) const
2037     {
2038     if (isExpanded()){
2039     right.expand();
2040     } else if (isTagged()) {
2041     if (right.isConstant()) {
2042     right.tag();
2043     }
2044     }
2045     }
2046    
2047     void
2048     Data::typeMatchRight(const Data& right)
2049     {
2050 jgs 94 if (isTagged()) {
2051     if (right.isExpanded()) {
2052     expand();
2053     }
2054     } else if (isConstant()) {
2055     if (right.isExpanded()) {
2056     expand();
2057     } else if (right.isTagged()) {
2058     tag();
2059     }
2060     }
2061     }
2062    
2063 bcumming 751 /* TODO */
2064     /* global reduction */
2065 jgs 94 void
2066     Data::setTaggedValue(int tagKey,
2067     const boost::python::object& value)
2068     {
2069 gross 783 if (isProtected()) {
2070     throw DataException("Error - attempt to update protected Data object.");
2071     }
2072 jgs 94 //
2073     // Ensure underlying data object is of type DataTagged
2074     tag();
2075    
2076     if (!isTagged()) {
2077     throw DataException("Error - DataTagged conversion failed!!");
2078     }
2079    
2080     //
2081     // Construct DataArray from boost::python::object input value
2082     DataArray valueDataArray(value);
2083    
2084     //
2085     // Call DataAbstract::setTaggedValue
2086     m_data->setTaggedValue(tagKey,valueDataArray.getView());
2087     }
2088    
2089 bcumming 751 /* TODO */
2090     /* global reduction */
2091 jgs 110 void
2092 jgs 121 Data::setTaggedValueFromCPP(int tagKey,
2093     const DataArrayView& value)
2094     {
2095 gross 783 if (isProtected()) {
2096     throw DataException("Error - attempt to update protected Data object.");
2097     }
2098 jgs 121 //
2099     // Ensure underlying data object is of type DataTagged
2100     tag();
2101    
2102     if (!isTagged()) {
2103     throw DataException("Error - DataTagged conversion failed!!");
2104     }
2105    
2106     //
2107     // Call DataAbstract::setTaggedValue
2108     m_data->setTaggedValue(tagKey,value);
2109     }
2110    
2111 bcumming 751 /* TODO */
2112     /* global reduction */
2113 jgs 149 int
2114     Data::getTagNumber(int dpno)
2115     {
2116     return m_data->getTagNumber(dpno);
2117     }
2118    
2119 bcumming 751 /* TODO */
2120     /* global reduction */
2121 jgs 121 void
2122 jgs 110 Data::setRefValue(int ref,
2123     const boost::python::numeric::array& value)
2124     {
2125 gross 783 if (isProtected()) {
2126     throw DataException("Error - attempt to update protected Data object.");
2127     }
2128 jgs 110 //
2129     // Construct DataArray from boost::python::object input value
2130     DataArray valueDataArray(value);
2131    
2132     //
2133     // Call DataAbstract::setRefValue
2134     m_data->setRefValue(ref,valueDataArray);
2135     }
2136    
2137 bcumming 751 /* TODO */
2138     /* global reduction */
2139 jgs 110 void
2140     Data::getRefValue(int ref,
2141     boost::python::numeric::array& value)
2142     {
2143     //
2144     // Construct DataArray for boost::python::object return value
2145     DataArray valueDataArray(value);
2146    
2147     //
2148     // Load DataArray with values from data-points specified by ref
2149     m_data->getRefValue(ref,valueDataArray);
2150    
2151     //
2152     // Load values from valueDataArray into return numarray
2153    
2154     // extract the shape of the numarray
2155     int rank = value.getrank();
2156     DataArrayView::ShapeType shape;
2157     for (int i=0; i < rank; i++) {
2158     shape.push_back(extract<int>(value.getshape()[i]));
2159     }
2160    
2161     // and load the numarray with the data from the DataArray
2162     DataArrayView valueView = valueDataArray.getView();
2163    
2164     if (rank==0) {
2165 jgs 126 boost::python::numeric::array temp_numArray(valueView());
2166     value = temp_numArray;
2167 jgs 110 }
2168     if (rank==1) {
2169     for (int i=0; i < shape[0]; i++) {
2170     value[i] = valueView(i);
2171     }
2172     }
2173     if (rank==2) {
2174 jgs 126 for (int i=0; i < shape[0]; i++) {
2175     for (int j=0; j < shape[1]; j++) {
2176     value[i][j] = valueView(i,j);
2177     }
2178     }
2179 jgs 110 }
2180     if (rank==3) {
2181 jgs 126 for (int i=0; i < shape[0]; i++) {
2182     for (int j=0; j < shape[1]; j++) {
2183     for (int k=0; k < shape[2]; k++) {
2184     value[i][j][k] = valueView(i,j,k);
2185     }
2186     }
2187     }
2188 jgs 110 }
2189     if (rank==4) {
2190 jgs 126 for (int i=0; i < shape[0]; i++) {
2191     for (int j=0; j < shape[1]; j++) {
2192     for (int k=0; k < shape[2]; k++) {
2193     for (int l=0; l < shape[3]; l++) {
2194     value[i][j][k][l] = valueView(i,j,k,l);
2195     }
2196     }
2197     }
2198     }
2199 jgs 110 }
2200    
2201     }
2202    
2203 jgs 94 void
2204 jgs 119 Data::archiveData(const std::string fileName)
2205     {
2206     cout << "Archiving Data object to: " << fileName << endl;
2207    
2208     //
2209     // Determine type of this Data object
2210     int dataType = -1;
2211    
2212     if (isEmpty()) {
2213     dataType = 0;
2214     cout << "\tdataType: DataEmpty" << endl;
2215     }
2216     if (isConstant()) {
2217     dataType = 1;
2218     cout << "\tdataType: DataConstant" << endl;
2219     }
2220     if (isTagged()) {
2221     dataType = 2;
2222     cout << "\tdataType: DataTagged" << endl;
2223     }
2224     if (isExpanded()) {
2225     dataType = 3;
2226     cout << "\tdataType: DataExpanded" << endl;
2227     }
2228 jgs 123
2229 jgs 119 if (dataType == -1) {
2230     throw DataException("archiveData Error: undefined dataType");
2231     }
2232    
2233     //
2234     // Collect data items common to all Data types
2235     int noSamples = getNumSamples();
2236     int noDPPSample = getNumDataPointsPerSample();
2237     int functionSpaceType = getFunctionSpace().getTypeCode();
2238     int dataPointRank = getDataPointRank();
2239     int dataPointSize = getDataPointSize();
2240     int dataLength = getLength();
2241     DataArrayView::ShapeType dataPointShape = getDataPointShape();
2242 woo409 757 vector<int> referenceNumbers(noSamples);
2243 jgs 119 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2244     referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
2245     }
2246 woo409 757 vector<int> tagNumbers(noSamples);
2247 jgs 119 if (isTagged()) {
2248     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2249     tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
2250     }
2251     }
2252    
2253     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2254     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2255     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2256    
2257     //
2258     // Flatten Shape to an array of integers suitable for writing to file
2259     int flatShape[4] = {0,0,0,0};
2260     cout << "\tshape: < ";
2261     for (int dim=0; dim<dataPointRank; dim++) {
2262     flatShape[dim] = dataPointShape[dim];
2263     cout << dataPointShape[dim] << " ";
2264     }
2265     cout << ">" << endl;
2266    
2267     //
2268 jgs 123 // Open archive file
2269 jgs 119 ofstream archiveFile;
2270     archiveFile.open(fileName.data(), ios::out);
2271    
2272     if (!archiveFile.good()) {
2273     throw DataException("archiveData Error: problem opening archive file");
2274     }
2275    
2276 jgs 123 //
2277     // Write common data items to archive file
2278 jgs 119 archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
2279     archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
2280     archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2281     archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2282     archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2283     archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2284     archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
2285     for (int dim = 0; dim < 4; dim++) {
2286     archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2287     }
2288     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2289     archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2290     }
2291     if (isTagged()) {
2292     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2293     archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2294     }
2295     }
2296    
2297     if (!archiveFile.good()) {
2298     throw DataException("archiveData Error: problem writing to archive file");
2299     }
2300    
2301     //
2302 jgs 123 // Archive underlying data values for each Data type
2303     int noValues;
2304 jgs 119 switch (dataType) {
2305     case 0:
2306     // DataEmpty
2307 jgs 123 noValues = 0;
2308     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2309     cout << "\tnoValues: " << noValues << endl;
2310 jgs 119 break;
2311     case 1:
2312     // DataConstant
2313 jgs 123 noValues = m_data->getLength();
2314     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2315     cout << "\tnoValues: " << noValues << endl;
2316     if (m_data->archiveData(archiveFile,noValues)) {
2317     throw DataException("archiveData Error: problem writing data to archive file");
2318     }
2319 jgs 119 break;
2320     case 2:
2321     // DataTagged
2322 jgs 123 noValues = m_data->getLength();
2323     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2324     cout << "\tnoValues: " << noValues << endl;
2325     if (m_data->archiveData(archiveFile,noValues)) {
2326     throw DataException("archiveData Error: problem writing data to archive file");
2327     }
2328 jgs 119 break;
2329     case 3:
2330     // DataExpanded
2331 jgs 123 noValues = m_data->getLength();
2332     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2333     cout << "\tnoValues: " << noValues << endl;
2334     if (m_data->archiveData(archiveFile,noValues)) {
2335     throw DataException("archiveData Error: problem writing data to archive file");
2336     }
2337 jgs 119 break;
2338     }
2339    
2340 jgs 123 if (!archiveFile.good()) {
2341     throw DataException("archiveData Error: problem writing data to archive file");
2342     }
2343    
2344     //
2345     // Close archive file
2346     archiveFile.close();
2347    
2348     if (!archiveFile.good()) {
2349     throw DataException("archiveData Error: problem closing archive file");
2350     }
2351    
2352 jgs 119 }
2353    
2354     void
2355     Data::extractData(const std::string fileName,
2356     const FunctionSpace& fspace)
2357     {
2358     //
2359     // Can only extract Data to an object which is initially DataEmpty
2360     if (!isEmpty()) {
2361     throw DataException("extractData Error: can only extract to DataEmpty object");
2362     }
2363    
2364     cout << "Extracting Data object from: " << fileName << endl;
2365    
2366     int dataType;
2367     int noSamples;
2368     int noDPPSample;
2369     int functionSpaceType;
2370     int dataPointRank;
2371     int dataPointSize;
2372     int dataLength;
2373     DataArrayView::ShapeType dataPointShape;
2374     int flatShape[4];
2375    
2376     //
2377 jgs 123 // Open the archive file
2378 jgs 119 ifstream archiveFile;
2379     archiveFile.open(fileName.data(), ios::in);
2380    
2381     if (!archiveFile.good()) {
2382     throw DataException("extractData Error: problem opening archive file");
2383     }
2384    
2385 jgs 123 //
2386     // Read common data items from archive file
2387 jgs 119 archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
2388     archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2389     archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2390     archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2391     archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2392     archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2393     archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
2394     for (int dim = 0; dim < 4; dim++) {
2395     archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2396     if (flatShape[dim]>0) {
2397     dataPointShape.push_back(flatShape[dim]);
2398     }
2399     }
2400 woo409 757 vector<int> referenceNumbers(noSamples);
2401 jgs 119 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2402     archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2403     }
2404 woo409 757 vector<int> tagNumbers(noSamples);
2405 jgs 119 if (dataType==2) {
2406     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2407     archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2408     }
2409     }
2410    
2411     if (!archiveFile.good()) {
2412     throw DataException("extractData Error: problem reading from archive file");
2413     }
2414    
2415 jgs 123 //
2416     // Verify the values just read from the archive file
2417 jgs 119 switch (dataType) {
2418     case 0:
2419     cout << "\tdataType: DataEmpty" << endl;
2420     break;
2421     case 1:
2422     cout << "\tdataType: DataConstant" << endl;
2423     break;
2424     case 2:
2425     cout << "\tdataType: DataTagged" << endl;
2426     break;
2427     case 3:
2428     cout << "\tdataType: DataExpanded" << endl;
2429     break;
2430     default:
2431     throw DataException("extractData Error: undefined dataType read from archive file");
2432     break;
2433     }
2434    
2435     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2436     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2437     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2438     cout << "\tshape: < ";
2439     for (int dim = 0; dim < dataPointRank; dim++) {
2440     cout << dataPointShape[dim] << " ";
2441     }
2442     cout << ">" << endl;
2443    
2444     //
2445     // Verify that supplied FunctionSpace object is compatible with this Data object.
2446     if ( (fspace.getTypeCode()!=functionSpaceType) ||
2447     (fspace.getNumSamples()!=noSamples) ||
2448     (fspace.getNumDPPSample()!=noDPPSample)
2449     ) {
2450     throw DataException("extractData Error: incompatible FunctionSpace");
2451     }
2452     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2453     if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2454     throw DataException("extractData Error: incompatible FunctionSpace");
2455     }
2456     }
2457     if (dataType==2) {
2458     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2459     if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2460     throw DataException("extractData Error: incompatible FunctionSpace");
2461     }
2462     }
2463     }
2464    
2465     //
2466     // Construct a DataVector to hold underlying data values
2467     DataVector dataVec(dataLength);
2468    
2469     //
2470     // Load this DataVector with the appropriate values
2471 jgs 123 int noValues;
2472     archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2473     cout << "\tnoValues: " << noValues << endl;
2474 jgs 119 switch (dataType) {
2475     case 0:
2476     // DataEmpty
2477 jgs 123 if (noValues != 0) {
2478     throw DataException("extractData Error: problem reading data from archive file");
2479     }
2480 jgs 119 break;
2481     case 1:
2482     // DataConstant
2483 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2484     throw DataException("extractData Error: problem reading data from archive file");
2485     }
2486 jgs 119 break;
2487     case 2:
2488     // DataTagged
2489 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2490     throw DataException("extractData Error: problem reading data from archive file");
2491     }
2492 jgs 119 break;
2493     case 3:
2494     // DataExpanded
2495 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2496     throw DataException("extractData Error: problem reading data from archive file");
2497     }
2498 jgs 119 break;
2499     }
2500    
2501 jgs 123 if (!archiveFile.good()) {
2502     throw DataException("extractData Error: problem reading from archive file");
2503     }
2504    
2505 jgs 119 //
2506 jgs 123 // Close archive file
2507     archiveFile.close();
2508    
2509     if (!archiveFile.good()) {
2510     throw DataException("extractData Error: problem closing archive file");
2511     }
2512    
2513     //
2514 jgs 119 // Construct an appropriate Data object
2515     DataAbstract* tempData;
2516     switch (dataType) {
2517     case 0:
2518     // DataEmpty
2519     tempData=new DataEmpty();
2520     break;
2521     case 1:
2522     // DataConstant
2523     tempData=new DataConstant(fspace,dataPointShape,dataVec);
2524     break;
2525     case 2:
2526     // DataTagged
2527     tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2528     break;
2529     case 3:
2530     // DataExpanded
2531     tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2532     break;
2533     }
2534     shared_ptr<DataAbstract> temp_data(tempData);
2535     m_data=temp_data;
2536     }
2537    
2538 jgs 94 ostream& escript::operator<<(ostream& o, const Data& data)
2539     {
2540     o << data.toString();
2541     return o;
2542     }
2543 bcumming 782
2544     /* Member functions specific to the MPI implementation */
2545