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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 921 - (hide annotations)
Fri Jan 5 00:54:37 2007 UTC (12 years, 8 months ago) by gross
File size: 83284 byte(s)
I have done some clarification on functions that allow to access individual data point values in a Data object. 
The term "data point number" is always local on a MPI process and referes to the value (data_point_in_sample, sample)
as a single identifyer (data_point_in_sample + sample * number_data_points_per_sample). a "global data point number"
referes to a tuple of a processour id and local data point number.

The function convertToNumArrayFromSampleNo has been removed now and convertToNumArrayFromDPNo renamed to getValueOfDataPoint.
There are two new functions:

   getNumberOfDataPoints
   setValueOfDataPoint

This allows you to do things like:

  in=Data(..)
  out=Data(..)
   for i in xrange(in.getNumberOfDataPoints())
       in_loc=in.getValueOfDataPoint(i)
       out_loc=< some operations on in_loc>
       out.setValueOfDataPoint(i,out_loc)


Also mindp  is renamed to  minGlobalDataPoint and there is a new function getValueOfGlobalDataPoint. While in MPI the functions getNumberOfDataPoints and getValueOfDataPoint are working locally on each process (so the code above is executed in parallel).
the latter allows getting a single value across all processors. 


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