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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 813 - (hide annotations)
Mon Aug 21 02:08:47 2006 UTC (13 years ago) by ksteube
File size: 83823 byte(s)
Tensor products for Data objects are now computed by a C++ method
C_GeneralTensorProduct, which calls C function matrix_matrix_product
to do the actual calculation.

Can perform product with either input transposed in place, meaning
without first computing the transpose in a separate step.

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 gross 804 Data::swapaxes(const int axis0, const int axis1) const
1320 jgs 123 {
1321 gross 804 int axis0_tmp,axis1_tmp;
1322 gross 800 #if defined DOPROF
1323     profData->unary++;
1324     #endif
1325     DataArrayView::ShapeType s=getDataPointShape();
1326     DataArrayView::ShapeType ev_shape;
1327     // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1328     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1329     int rank=getDataPointRank();
1330 gross 804 if (rank<2) {
1331     throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1332 gross 800 }
1333 gross 804 if (axis0<0 || axis0>rank-1) {
1334     throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1335     }
1336     if (axis1<0 || axis1>rank-1) {
1337     throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1338     }
1339     if (axis0 == axis1) {
1340     throw DataException("Error - Data::swapaxes: axis indices must be different.");
1341     }
1342     if (axis0 > axis1) {
1343     axis0_tmp=axis1;
1344     axis1_tmp=axis0;
1345     } else {
1346     axis0_tmp=axis0;
1347     axis1_tmp=axis1;
1348     }
1349 gross 800 for (int i=0; i<rank; i++) {
1350 gross 804 if (i == axis0_tmp) {
1351     ev_shape.push_back(s[axis1_tmp]);
1352     } else if (i == axis1_tmp) {
1353     ev_shape.push_back(s[axis0_tmp]);
1354 gross 800 } else {
1355     ev_shape.push_back(s[i]);
1356     }
1357     }
1358     Data ev(0.,ev_shape,getFunctionSpace());
1359     ev.typeMatchRight(*this);
1360 gross 804 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1361 gross 800 return ev;
1362    
1363 jgs 123 }
1364    
1365     Data
1366 ksteube 775 Data::symmetric() const
1367 jgs 123 {
1368 ksteube 775 #if defined DOPROF
1369     profData->unary++;
1370     #endif
1371     // check input
1372     DataArrayView::ShapeType s=getDataPointShape();
1373     if (getDataPointRank()==2) {
1374     if(s[0] != s[1])
1375     throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1376     }
1377     else if (getDataPointRank()==4) {
1378     if(!(s[0] == s[2] && s[1] == s[3]))
1379     throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1380     }
1381     else {
1382     throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1383     }
1384     Data ev(0.,getDataPointShape(),getFunctionSpace());
1385     ev.typeMatchRight(*this);
1386     m_data->symmetric(ev.m_data.get());
1387     return ev;
1388     }
1389    
1390     Data
1391     Data::nonsymmetric() const
1392     {
1393     #if defined DOPROF
1394     profData->unary++;
1395     #endif
1396     // check input
1397     DataArrayView::ShapeType s=getDataPointShape();
1398     if (getDataPointRank()==2) {
1399     if(s[0] != s[1])
1400     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1401     DataArrayView::ShapeType ev_shape;
1402     ev_shape.push_back(s[0]);
1403     ev_shape.push_back(s[1]);
1404     Data ev(0.,ev_shape,getFunctionSpace());
1405     ev.typeMatchRight(*this);
1406     m_data->nonsymmetric(ev.m_data.get());
1407     return ev;
1408     }
1409     else if (getDataPointRank()==4) {
1410     if(!(s[0] == s[2] && s[1] == s[3]))
1411     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1412     DataArrayView::ShapeType ev_shape;
1413     ev_shape.push_back(s[0]);
1414     ev_shape.push_back(s[1]);
1415     ev_shape.push_back(s[2]);
1416     ev_shape.push_back(s[3]);
1417     Data ev(0.,ev_shape,getFunctionSpace());
1418     ev.typeMatchRight(*this);
1419     m_data->nonsymmetric(ev.m_data.get());
1420     return ev;
1421     }
1422     else {
1423     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1424     }
1425     }
1426    
1427     Data
1428 gross 800 Data::trace(int axis_offset) const
1429 ksteube 775 {
1430     #if defined DOPROF
1431     profData->unary++;
1432     #endif
1433     DataArrayView::ShapeType s=getDataPointShape();
1434     if (getDataPointRank()==2) {
1435     DataArrayView::ShapeType ev_shape;
1436     Data ev(0.,ev_shape,getFunctionSpace());
1437     ev.typeMatchRight(*this);
1438 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1439 ksteube 775 return ev;
1440     }
1441     if (getDataPointRank()==3) {
1442     DataArrayView::ShapeType ev_shape;
1443     if (axis_offset==0) {
1444     int s2=s[2];
1445     ev_shape.push_back(s2);
1446     }
1447     else if (axis_offset==1) {
1448     int s0=s[0];
1449     ev_shape.push_back(s0);
1450     }
1451     Data ev(0.,ev_shape,getFunctionSpace());
1452     ev.typeMatchRight(*this);
1453 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1454 ksteube 775 return ev;
1455     }
1456     if (getDataPointRank()==4) {
1457     DataArrayView::ShapeType ev_shape;
1458     if (axis_offset==0) {
1459     ev_shape.push_back(s[2]);
1460     ev_shape.push_back(s[3]);
1461     }
1462     else if (axis_offset==1) {
1463     ev_shape.push_back(s[0]);
1464     ev_shape.push_back(s[3]);
1465     }
1466     else if (axis_offset==2) {
1467     ev_shape.push_back(s[0]);
1468     ev_shape.push_back(s[1]);
1469     }
1470     Data ev(0.,ev_shape,getFunctionSpace());
1471     ev.typeMatchRight(*this);
1472 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1473 ksteube 775 return ev;
1474     }
1475     else {
1476 gross 800 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1477 ksteube 775 }
1478     }
1479    
1480     Data
1481     Data::transpose(int axis_offset) const
1482     {
1483 gross 800 #if defined DOPROF
1484     profData->unary++;
1485     #endif
1486 ksteube 775 DataArrayView::ShapeType s=getDataPointShape();
1487     DataArrayView::ShapeType ev_shape;
1488     // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1489     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1490     int rank=getDataPointRank();
1491     if (axis_offset<0 || axis_offset>rank) {
1492     throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1493     }
1494     for (int i=0; i<rank; i++) {
1495     int index = (axis_offset+i)%rank;
1496     ev_shape.push_back(s[index]); // Append to new shape
1497     }
1498     Data ev(0.,ev_shape,getFunctionSpace());
1499     ev.typeMatchRight(*this);
1500     m_data->transpose(ev.m_data.get(), axis_offset);
1501     return ev;
1502 jgs 123 }
1503    
1504 gross 576 Data
1505     Data::eigenvalues() const
1506     {
1507     #if defined DOPROF
1508     profData->unary++;
1509     #endif
1510     // check input
1511     DataArrayView::ShapeType s=getDataPointShape();
1512     if (getDataPointRank()!=2)
1513     throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1514     if(s[0] != s[1])
1515     throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1516     // create return
1517     DataArrayView::ShapeType ev_shape(1,s[0]);
1518     Data ev(0.,ev_shape,getFunctionSpace());
1519     ev.typeMatchRight(*this);
1520     m_data->eigenvalues(ev.m_data.get());
1521     return ev;
1522     }
1523    
1524 jgs 121 const boost::python::tuple
1525 gross 576 Data::eigenvalues_and_eigenvectors(const double tol) const
1526     {
1527     #if defined DOPROF
1528     profData->unary++;
1529     #endif
1530     DataArrayView::ShapeType s=getDataPointShape();
1531     if (getDataPointRank()!=2)
1532     throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1533     if(s[0] != s[1])
1534     throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1535     // create return
1536     DataArrayView::ShapeType ev_shape(1,s[0]);
1537     Data ev(0.,ev_shape,getFunctionSpace());
1538     ev.typeMatchRight(*this);
1539     DataArrayView::ShapeType V_shape(2,s[0]);
1540     Data V(0.,V_shape,getFunctionSpace());
1541     V.typeMatchRight(*this);
1542     m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1543     return make_tuple(boost::python::object(ev),boost::python::object(V));
1544     }
1545    
1546     const boost::python::tuple
1547 jgs 121 Data::mindp() const
1548     {
1549 jgs 148 // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1550     // abort (for unknown reasons) if there are openmp directives with it in the
1551     // surrounding function
1552    
1553     int SampleNo;
1554     int DataPointNo;
1555 bcumming 782 int ProcNo;
1556     calc_mindp(ProcNo,SampleNo,DataPointNo);
1557     return make_tuple(ProcNo,SampleNo,DataPointNo);
1558 jgs 148 }
1559    
1560     void
1561 bcumming 782 Data::calc_mindp( int& ProcNo,
1562     int& SampleNo,
1563     int& DataPointNo) const
1564 jgs 148 {
1565     int i,j;
1566     int lowi=0,lowj=0;
1567     double min=numeric_limits<double>::max();
1568    
1569 jgs 121 Data temp=minval();
1570    
1571     int numSamples=temp.getNumSamples();
1572     int numDPPSample=temp.getNumDataPointsPerSample();
1573    
1574 jgs 148 double next,local_min;
1575     int local_lowi,local_lowj;
1576 jgs 121
1577 jgs 148 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1578     {
1579     local_min=min;
1580     #pragma omp for private(i,j) schedule(static)
1581     for (i=0; i<numSamples; i++) {
1582     for (j=0; j<numDPPSample; j++) {
1583     next=temp.getDataPoint(i,j)();
1584     if (next<local_min) {
1585     local_min=next;
1586     local_lowi=i;
1587     local_lowj=j;
1588     }
1589 jgs 121 }
1590     }
1591 jgs 148 #pragma omp critical
1592     if (local_min<min) {
1593     min=local_min;
1594     lowi=local_lowi;
1595     lowj=local_lowj;
1596     }
1597 jgs 121 }
1598    
1599 bcumming 782 #ifdef PASO_MPI
1600     // determine the processor on which the minimum occurs
1601     next = temp.getDataPoint(lowi,lowj)();
1602     int lowProc = 0;
1603     double *globalMins = new double[get_MPISize()+1];
1604     int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1605    
1606     if( get_MPIRank()==0 ){
1607     next = globalMins[lowProc];
1608     for( i=1; i<get_MPISize(); i++ )
1609     if( next>globalMins[i] ){
1610     lowProc = i;
1611     next = globalMins[i];
1612     }
1613     }
1614     MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1615    
1616     delete [] globalMins;
1617     ProcNo = lowProc;
1618 bcumming 790 #else
1619     ProcNo = 0;
1620 bcumming 782 #endif
1621 jgs 148 SampleNo = lowi;
1622     DataPointNo = lowj;
1623 jgs 121 }
1624    
1625 jgs 104 void
1626     Data::saveDX(std::string fileName) const
1627     {
1628 jgs 153 boost::python::dict args;
1629     args["data"]=boost::python::object(this);
1630     getDomain().saveDX(fileName,args);
1631 jgs 104 return;
1632     }
1633    
1634 jgs 110 void
1635     Data::saveVTK(std::string fileName) const
1636     {
1637 jgs 153 boost::python::dict args;
1638     args["data"]=boost::python::object(this);
1639     getDomain().saveVTK(fileName,args);
1640 jgs 110 return;
1641     }
1642    
1643 jgs 102 Data&
1644     Data::operator+=(const Data& right)
1645     {
1646 gross 783 if (isProtected()) {
1647     throw DataException("Error - attempt to update protected Data object.");
1648     }
1649 jgs 151 #if defined DOPROF
1650 jgs 123 profData->binary++;
1651 jgs 151 #endif
1652 jgs 94 binaryOp(right,plus<double>());
1653     return (*this);
1654     }
1655    
1656 jgs 102 Data&
1657     Data::operator+=(const boost::python::object& right)
1658 jgs 94 {
1659 gross 783 if (isProtected()) {
1660     throw DataException("Error - attempt to update protected Data object.");
1661     }
1662 jgs 151 #if defined DOPROF
1663 jgs 123 profData->binary++;
1664 jgs 151 #endif
1665 jgs 94 binaryOp(right,plus<double>());
1666     return (*this);
1667     }
1668    
1669 jgs 102 Data&
1670     Data::operator-=(const Data& right)
1671 jgs 94 {
1672 gross 783 if (isProtected()) {
1673     throw DataException("Error - attempt to update protected Data object.");
1674     }
1675 jgs 151 #if defined DOPROF
1676 jgs 123 profData->binary++;
1677 jgs 151 #endif
1678 jgs 94 binaryOp(right,minus<double>());
1679     return (*this);
1680     }
1681    
1682 jgs 102 Data&
1683     Data::operator-=(const boost::python::object& right)
1684 jgs 94 {
1685 gross 783 if (isProtected()) {
1686     throw DataException("Error - attempt to update protected Data object.");
1687     }
1688 jgs 151 #if defined DOPROF
1689 jgs 123 profData->binary++;
1690 jgs 151 #endif
1691 jgs 94 binaryOp(right,minus<double>());
1692     return (*this);
1693     }
1694    
1695 jgs 102 Data&
1696     Data::operator*=(const Data& right)
1697 jgs 94 {
1698 gross 783 if (isProtected()) {
1699     throw DataException("Error - attempt to update protected Data object.");
1700     }
1701 jgs 151 #if defined DOPROF
1702 jgs 123 profData->binary++;
1703 jgs 151 #endif
1704 jgs 94 binaryOp(right,multiplies<double>());
1705     return (*this);
1706     }
1707    
1708 jgs 102 Data&
1709     Data::operator*=(const boost::python::object& right)
1710 jgs 94 {
1711 gross 783 if (isProtected()) {
1712     throw DataException("Error - attempt to update protected Data object.");
1713     }
1714 jgs 151 #if defined DOPROF
1715 jgs 123 profData->binary++;
1716 jgs 151 #endif
1717 jgs 94 binaryOp(right,multiplies<double>());
1718     return (*this);
1719     }
1720    
1721 jgs 102 Data&
1722     Data::operator/=(const Data& right)
1723 jgs 94 {
1724 gross 783 if (isProtected()) {
1725     throw DataException("Error - attempt to update protected Data object.");
1726     }
1727 jgs 151 #if defined DOPROF
1728 jgs 123 profData->binary++;
1729 jgs 151 #endif
1730 jgs 94 binaryOp(right,divides<double>());
1731     return (*this);
1732     }
1733    
1734 jgs 102 Data&
1735     Data::operator/=(const boost::python::object& right)
1736 jgs 94 {
1737 gross 783 if (isProtected()) {
1738     throw DataException("Error - attempt to update protected Data object.");
1739     }
1740 jgs 151 #if defined DOPROF
1741 jgs 123 profData->binary++;
1742 jgs 151 #endif
1743 jgs 94 binaryOp(right,divides<double>());
1744     return (*this);
1745     }
1746    
1747 jgs 102 Data
1748 gross 699 Data::rpowO(const boost::python::object& left) const
1749     {
1750 gross 783 if (isProtected()) {
1751     throw DataException("Error - attempt to update protected Data object.");
1752     }
1753 gross 699 #if defined DOPROF
1754     profData->binary++;
1755     #endif
1756     Data left_d(left,*this);
1757     return left_d.powD(*this);
1758     }
1759    
1760     Data
1761 jgs 102 Data::powO(const boost::python::object& right) const
1762 jgs 94 {
1763 jgs 151 #if defined DOPROF
1764 jgs 123 profData->binary++;
1765 jgs 151 #endif
1766 jgs 94 Data result;
1767     result.copy(*this);
1768     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1769     return result;
1770     }
1771    
1772 jgs 102 Data
1773     Data::powD(const Data& right) const
1774 jgs 94 {
1775 jgs 151 #if defined DOPROF
1776 jgs 123 profData->binary++;
1777 jgs 151 #endif
1778 jgs 94 Data result;
1779     result.copy(*this);
1780     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1781     return result;
1782     }
1783    
1784 bcumming 751
1785 jgs 94 //
1786 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1787 jgs 102 Data
1788     escript::operator+(const Data& left, const Data& right)
1789 jgs 94 {
1790     Data result;
1791     //
1792     // perform a deep copy
1793     result.copy(left);
1794     result+=right;
1795     return result;
1796     }
1797    
1798     //
1799 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1800 jgs 102 Data
1801     escript::operator-(const Data& left, const Data& right)
1802 jgs 94 {
1803     Data result;
1804     //
1805     // perform a deep copy
1806     result.copy(left);
1807     result-=right;
1808     return result;
1809     }
1810    
1811     //
1812 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1813 jgs 102 Data
1814     escript::operator*(const Data& left, const Data& right)
1815 jgs 94 {
1816     Data result;
1817     //
1818     // perform a deep copy
1819     result.copy(left);
1820     result*=right;
1821     return result;
1822     }
1823    
1824     //
1825 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1826 jgs 102 Data
1827     escript::operator/(const Data& left, const Data& right)
1828 jgs 94 {
1829     Data result;
1830     //
1831     // perform a deep copy
1832     result.copy(left);
1833     result/=right;
1834     return result;
1835     }
1836    
1837     //
1838 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1839 jgs 102 Data
1840     escript::operator+(const Data& left, const boost::python::object& right)
1841 jgs 94 {
1842     //
1843     // Convert to DataArray format if possible
1844     DataArray temp(right);
1845     Data result;
1846     //
1847     // perform a deep copy
1848     result.copy(left);
1849     result+=right;
1850     return result;
1851     }
1852    
1853     //
1854 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1855 jgs 102 Data
1856     escript::operator-(const Data& left, const boost::python::object& right)
1857 jgs 94 {
1858     //
1859     // Convert to DataArray format if possible
1860     DataArray temp(right);
1861     Data result;
1862     //
1863     // perform a deep copy
1864     result.copy(left);
1865     result-=right;
1866     return result;
1867     }
1868    
1869     //
1870 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1871 jgs 102 Data
1872     escript::operator*(const Data& left, const boost::python::object& right)
1873 jgs 94 {
1874     //
1875     // Convert to DataArray format if possible
1876     DataArray temp(right);
1877     Data result;
1878     //
1879     // perform a deep copy
1880     result.copy(left);
1881     result*=right;
1882     return result;
1883     }
1884    
1885     //
1886 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1887 jgs 102 Data
1888     escript::operator/(const Data& left, const boost::python::object& right)
1889 jgs 94 {
1890     //
1891     // Convert to DataArray format if possible
1892     DataArray temp(right);
1893     Data result;
1894     //
1895     // perform a deep copy
1896     result.copy(left);
1897     result/=right;
1898     return result;
1899     }
1900    
1901     //
1902 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1903 jgs 102 Data
1904     escript::operator+(const boost::python::object& left, const Data& right)
1905 jgs 94 {
1906     //
1907     // Construct the result using the given value and the other parameters
1908     // from right
1909     Data result(left,right);
1910     result+=right;
1911     return result;
1912     }
1913    
1914     //
1915 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1916 jgs 102 Data
1917     escript::operator-(const boost::python::object& left, const Data& right)
1918 jgs 94 {
1919     //
1920     // Construct the result using the given value and the other parameters
1921     // from right
1922     Data result(left,right);
1923     result-=right;
1924     return result;
1925     }
1926    
1927     //
1928 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1929 jgs 102 Data
1930     escript::operator*(const boost::python::object& left, const Data& right)
1931 jgs 94 {
1932     //
1933     // Construct the result using the given value and the other parameters
1934     // from right
1935     Data result(left,right);
1936     result*=right;
1937     return result;
1938     }
1939    
1940     //
1941 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1942 jgs 102 Data
1943     escript::operator/(const boost::python::object& left, const Data& right)
1944 jgs 94 {
1945     //
1946     // Construct the result using the given value and the other parameters
1947     // from right
1948     Data result(left,right);
1949     result/=right;
1950     return result;
1951     }
1952    
1953     //
1954 jgs 102 //bool escript::operator==(const Data& left, const Data& right)
1955     //{
1956     // /*
1957     // NB: this operator does very little at this point, and isn't to
1958     // be relied on. Requires further implementation.
1959     // */
1960     //
1961     // bool ret;
1962     //
1963     // if (left.isEmpty()) {
1964     // if(!right.isEmpty()) {
1965     // ret = false;
1966     // } else {
1967     // ret = true;
1968     // }
1969     // }
1970     //
1971     // if (left.isConstant()) {
1972     // if(!right.isConstant()) {
1973     // ret = false;
1974     // } else {
1975     // ret = true;
1976     // }
1977     // }
1978     //
1979     // if (left.isTagged()) {
1980     // if(!right.isTagged()) {
1981     // ret = false;
1982     // } else {
1983     // ret = true;
1984     // }
1985     // }
1986     //
1987     // if (left.isExpanded()) {
1988     // if(!right.isExpanded()) {
1989     // ret = false;
1990     // } else {
1991     // ret = true;
1992     // }
1993     // }
1994     //
1995     // return ret;
1996     //}
1997    
1998 bcumming 751 /* TODO */
1999     /* global reduction */
2000 jgs 102 Data
2001     Data::getItem(const boost::python::object& key) const
2002 jgs 94 {
2003 jgs 102 const DataArrayView& view=getPointDataView();
2004 jgs 94
2005 jgs 102 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
2006 jgs 94
2007 jgs 102 if (slice_region.size()!=view.getRank()) {
2008     throw DataException("Error - slice size does not match Data rank.");
2009 jgs 94 }
2010    
2011 jgs 102 return getSlice(slice_region);
2012 jgs 94 }
2013    
2014 bcumming 751 /* TODO */
2015     /* global reduction */
2016 jgs 94 Data
2017 jgs 102 Data::getSlice(const DataArrayView::RegionType& region) const
2018 jgs 94 {
2019 jgs 151 #if defined DOPROF
2020 jgs 123 profData->slicing++;
2021 jgs 151 #endif
2022 jgs 102 return Data(*this,region);
2023 jgs 94 }
2024    
2025 bcumming 751 /* TODO */
2026     /* global reduction */
2027 jgs 94 void
2028 jgs 102 Data::setItemO(const boost::python::object& key,
2029     const boost::python::object& value)
2030 jgs 94 {
2031 jgs 102 Data tempData(value,getFunctionSpace());
2032     setItemD(key,tempData);
2033     }
2034    
2035 bcumming 751 /* TODO */
2036     /* global reduction */
2037 jgs 102 void
2038     Data::setItemD(const boost::python::object& key,
2039     const Data& value)
2040     {
2041 jgs 94 const DataArrayView& view=getPointDataView();
2042 jgs 123
2043 jgs 94 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
2044     if (slice_region.size()!=view.getRank()) {
2045     throw DataException("Error - slice size does not match Data rank.");
2046     }
2047 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
2048     setSlice(Data(value,getFunctionSpace()),slice_region);
2049     } else {
2050     setSlice(value,slice_region);
2051     }
2052 jgs 94 }
2053    
2054 bcumming 751 /* TODO */
2055     /* global reduction */
2056 jgs 94 void
2057 jgs 102 Data::setSlice(const Data& value,
2058     const DataArrayView::RegionType& region)
2059 jgs 94 {
2060 gross 783 if (isProtected()) {
2061     throw DataException("Error - attempt to update protected Data object.");
2062     }
2063 jgs 151 #if defined DOPROF
2064 jgs 123 profData->slicing++;
2065 jgs 151 #endif
2066 jgs 102 Data tempValue(value);
2067     typeMatchLeft(tempValue);
2068     typeMatchRight(tempValue);
2069     m_data->setSlice(tempValue.m_data.get(),region);
2070     }
2071    
2072     void
2073     Data::typeMatchLeft(Data& right) const
2074     {
2075     if (isExpanded()){
2076     right.expand();
2077     } else if (isTagged()) {
2078     if (right.isConstant()) {
2079     right.tag();
2080     }
2081     }
2082     }
2083    
2084     void
2085     Data::typeMatchRight(const Data& right)
2086     {
2087 jgs 94 if (isTagged()) {
2088     if (right.isExpanded()) {
2089     expand();
2090     }
2091     } else if (isConstant()) {
2092     if (right.isExpanded()) {
2093     expand();
2094     } else if (right.isTagged()) {
2095     tag();
2096     }
2097     }
2098     }
2099    
2100 bcumming 751 /* TODO */
2101     /* global reduction */
2102 jgs 94 void
2103     Data::setTaggedValue(int tagKey,
2104     const boost::python::object& value)
2105     {
2106 gross 783 if (isProtected()) {
2107     throw DataException("Error - attempt to update protected Data object.");
2108     }
2109 jgs 94 //
2110     // Ensure underlying data object is of type DataTagged
2111     tag();
2112    
2113     if (!isTagged()) {
2114     throw DataException("Error - DataTagged conversion failed!!");
2115     }
2116    
2117     //
2118     // Construct DataArray from boost::python::object input value
2119     DataArray valueDataArray(value);
2120    
2121     //
2122     // Call DataAbstract::setTaggedValue
2123     m_data->setTaggedValue(tagKey,valueDataArray.getView());
2124     }
2125    
2126 bcumming 751 /* TODO */
2127     /* global reduction */
2128 jgs 110 void
2129 jgs 121 Data::setTaggedValueFromCPP(int tagKey,
2130     const DataArrayView& value)
2131     {
2132 gross 783 if (isProtected()) {
2133     throw DataException("Error - attempt to update protected Data object.");
2134     }
2135 jgs 121 //
2136     // Ensure underlying data object is of type DataTagged
2137     tag();
2138    
2139     if (!isTagged()) {
2140     throw DataException("Error - DataTagged conversion failed!!");
2141     }
2142    
2143     //
2144     // Call DataAbstract::setTaggedValue
2145     m_data->setTaggedValue(tagKey,value);
2146     }
2147    
2148 bcumming 751 /* TODO */
2149     /* global reduction */
2150 jgs 149 int
2151     Data::getTagNumber(int dpno)
2152     {
2153     return m_data->getTagNumber(dpno);
2154     }
2155    
2156 bcumming 751 /* TODO */
2157     /* global reduction */
2158 jgs 121 void
2159 jgs 110 Data::setRefValue(int ref,
2160     const boost::python::numeric::array& value)
2161     {
2162 gross 783 if (isProtected()) {
2163     throw DataException("Error - attempt to update protected Data object.");
2164     }
2165 jgs 110 //
2166     // Construct DataArray from boost::python::object input value
2167     DataArray valueDataArray(value);
2168    
2169     //
2170     // Call DataAbstract::setRefValue
2171     m_data->setRefValue(ref,valueDataArray);
2172     }
2173    
2174 bcumming 751 /* TODO */
2175     /* global reduction */
2176 jgs 110 void
2177     Data::getRefValue(int ref,
2178     boost::python::numeric::array& value)
2179     {
2180     //
2181     // Construct DataArray for boost::python::object return value
2182     DataArray valueDataArray(value);
2183    
2184     //
2185     // Load DataArray with values from data-points specified by ref
2186     m_data->getRefValue(ref,valueDataArray);
2187    
2188     //
2189     // Load values from valueDataArray into return numarray
2190    
2191     // extract the shape of the numarray
2192     int rank = value.getrank();
2193     DataArrayView::ShapeType shape;
2194     for (int i=0; i < rank; i++) {
2195     shape.push_back(extract<int>(value.getshape()[i]));
2196     }
2197    
2198     // and load the numarray with the data from the DataArray
2199     DataArrayView valueView = valueDataArray.getView();
2200    
2201     if (rank==0) {
2202 jgs 126 boost::python::numeric::array temp_numArray(valueView());
2203     value = temp_numArray;
2204 jgs 110 }
2205     if (rank==1) {
2206     for (int i=0; i < shape[0]; i++) {
2207     value[i] = valueView(i);
2208     }
2209     }
2210     if (rank==2) {
2211 jgs 126 for (int i=0; i < shape[0]; i++) {
2212     for (int j=0; j < shape[1]; j++) {
2213     value[i][j] = valueView(i,j);
2214     }
2215     }
2216 jgs 110 }
2217     if (rank==3) {
2218 jgs 126 for (int i=0; i < shape[0]; i++) {
2219     for (int j=0; j < shape[1]; j++) {
2220     for (int k=0; k < shape[2]; k++) {
2221     value[i][j][k] = valueView(i,j,k);
2222     }
2223     }
2224     }
2225 jgs 110 }
2226     if (rank==4) {
2227 jgs 126 for (int i=0; i < shape[0]; i++) {
2228     for (int j=0; j < shape[1]; j++) {
2229     for (int k=0; k < shape[2]; k++) {
2230     for (int l=0; l < shape[3]; l++) {
2231     value[i][j][k][l] = valueView(i,j,k,l);
2232     }
2233     }
2234     }
2235     }
2236 jgs 110 }
2237    
2238     }
2239    
2240 jgs 94 void
2241 jgs 119 Data::archiveData(const std::string fileName)
2242     {
2243     cout << "Archiving Data object to: " << fileName << endl;
2244    
2245     //
2246     // Determine type of this Data object
2247     int dataType = -1;
2248    
2249     if (isEmpty()) {
2250     dataType = 0;
2251     cout << "\tdataType: DataEmpty" << endl;
2252     }
2253     if (isConstant()) {
2254     dataType = 1;
2255     cout << "\tdataType: DataConstant" << endl;
2256     }
2257     if (isTagged()) {
2258     dataType = 2;
2259     cout << "\tdataType: DataTagged" << endl;
2260     }
2261     if (isExpanded()) {
2262     dataType = 3;
2263     cout << "\tdataType: DataExpanded" << endl;
2264     }
2265 jgs 123
2266 jgs 119 if (dataType == -1) {
2267     throw DataException("archiveData Error: undefined dataType");
2268     }
2269    
2270     //
2271     // Collect data items common to all Data types
2272     int noSamples = getNumSamples();
2273     int noDPPSample = getNumDataPointsPerSample();
2274     int functionSpaceType = getFunctionSpace().getTypeCode();
2275     int dataPointRank = getDataPointRank();
2276     int dataPointSize = getDataPointSize();
2277     int dataLength = getLength();
2278     DataArrayView::ShapeType dataPointShape = getDataPointShape();
2279 woo409 757 vector<int> referenceNumbers(noSamples);
2280 jgs 119 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2281     referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
2282     }
2283 woo409 757 vector<int> tagNumbers(noSamples);
2284 jgs 119 if (isTagged()) {
2285     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2286     tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
2287     }
2288     }
2289    
2290     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2291     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2292     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2293    
2294     //
2295     // Flatten Shape to an array of integers suitable for writing to file
2296     int flatShape[4] = {0,0,0,0};
2297     cout << "\tshape: < ";
2298     for (int dim=0; dim<dataPointRank; dim++) {
2299     flatShape[dim] = dataPointShape[dim];
2300     cout << dataPointShape[dim] << " ";
2301     }
2302     cout << ">" << endl;
2303    
2304     //
2305 jgs 123 // Open archive file
2306 jgs 119 ofstream archiveFile;
2307     archiveFile.open(fileName.data(), ios::out);
2308    
2309     if (!archiveFile.good()) {
2310     throw DataException("archiveData Error: problem opening archive file");
2311     }
2312    
2313 jgs 123 //
2314     // Write common data items to archive file
2315 jgs 119 archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
2316     archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
2317     archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2318     archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2319     archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2320     archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2321     archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
2322     for (int dim = 0; dim < 4; dim++) {
2323     archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2324     }
2325     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2326     archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2327     }
2328     if (isTagged()) {
2329     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2330     archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2331     }
2332     }
2333    
2334     if (!archiveFile.good()) {
2335     throw DataException("archiveData Error: problem writing to archive file");
2336     }
2337    
2338     //
2339 jgs 123 // Archive underlying data values for each Data type
2340     int noValues;
2341 jgs 119 switch (dataType) {
2342     case 0:
2343     // DataEmpty
2344 jgs 123 noValues = 0;
2345     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2346     cout << "\tnoValues: " << noValues << endl;
2347 jgs 119 break;
2348     case 1:
2349     // DataConstant
2350 jgs 123 noValues = m_data->getLength();
2351     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2352     cout << "\tnoValues: " << noValues << endl;
2353     if (m_data->archiveData(archiveFile,noValues)) {
2354     throw DataException("archiveData Error: problem writing data to archive file");
2355     }
2356 jgs 119 break;
2357     case 2:
2358     // DataTagged
2359 jgs 123 noValues = m_data->getLength();
2360     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2361     cout << "\tnoValues: " << noValues << endl;
2362     if (m_data->archiveData(archiveFile,noValues)) {
2363     throw DataException("archiveData Error: problem writing data to archive file");
2364     }
2365 jgs 119 break;
2366     case 3:
2367     // DataExpanded
2368 jgs 123 noValues = m_data->getLength();
2369     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2370     cout << "\tnoValues: " << noValues << endl;
2371     if (m_data->archiveData(archiveFile,noValues)) {
2372     throw DataException("archiveData Error: problem writing data to archive file");
2373     }
2374 jgs 119 break;
2375     }
2376    
2377 jgs 123 if (!archiveFile.good()) {
2378     throw DataException("archiveData Error: problem writing data to archive file");
2379     }
2380    
2381     //
2382     // Close archive file
2383     archiveFile.close();
2384    
2385     if (!archiveFile.good()) {
2386     throw DataException("archiveData Error: problem closing archive file");
2387     }
2388    
2389 jgs 119 }
2390    
2391     void
2392     Data::extractData(const std::string fileName,
2393     const FunctionSpace& fspace)
2394     {
2395     //
2396     // Can only extract Data to an object which is initially DataEmpty
2397     if (!isEmpty()) {
2398     throw DataException("extractData Error: can only extract to DataEmpty object");
2399     }
2400    
2401     cout << "Extracting Data object from: " << fileName << endl;
2402    
2403     int dataType;
2404     int noSamples;
2405     int noDPPSample;
2406     int functionSpaceType;
2407     int dataPointRank;
2408     int dataPointSize;
2409     int dataLength;
2410     DataArrayView::ShapeType dataPointShape;
2411     int flatShape[4];
2412    
2413     //
2414 jgs 123 // Open the archive file
2415 jgs 119 ifstream archiveFile;
2416     archiveFile.open(fileName.data(), ios::in);
2417    
2418     if (!archiveFile.good()) {
2419     throw DataException("extractData Error: problem opening archive file");
2420     }
2421    
2422 jgs 123 //
2423     // Read common data items from archive file
2424 jgs 119 archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
2425     archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2426     archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2427     archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2428     archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2429     archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2430     archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
2431     for (int dim = 0; dim < 4; dim++) {
2432     archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2433     if (flatShape[dim]>0) {
2434     dataPointShape.push_back(flatShape[dim]);
2435     }
2436     }
2437 woo409 757 vector<int> referenceNumbers(noSamples);
2438 jgs 119 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2439     archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2440     }
2441 woo409 757 vector<int> tagNumbers(noSamples);
2442 jgs 119 if (dataType==2) {
2443     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2444     archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2445     }
2446     }
2447    
2448     if (!archiveFile.good()) {
2449     throw DataException("extractData Error: problem reading from archive file");
2450     }
2451    
2452 jgs 123 //
2453     // Verify the values just read from the archive file
2454 jgs 119 switch (dataType) {
2455     case 0:
2456     cout << "\tdataType: DataEmpty" << endl;
2457     break;
2458     case 1:
2459     cout << "\tdataType: DataConstant" << endl;
2460     break;
2461     case 2:
2462     cout << "\tdataType: DataTagged" << endl;
2463     break;
2464     case 3:
2465     cout << "\tdataType: DataExpanded" << endl;
2466     break;
2467     default:
2468     throw DataException("extractData Error: undefined dataType read from archive file");
2469     break;
2470     }
2471    
2472     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2473     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2474     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2475     cout << "\tshape: < ";
2476     for (int dim = 0; dim < dataPointRank; dim++) {
2477     cout << dataPointShape[dim] << " ";
2478     }
2479     cout << ">" << endl;
2480    
2481     //
2482     // Verify that supplied FunctionSpace object is compatible with this Data object.
2483     if ( (fspace.getTypeCode()!=functionSpaceType) ||
2484     (fspace.getNumSamples()!=noSamples) ||
2485     (fspace.getNumDPPSample()!=noDPPSample)
2486     ) {
2487     throw DataException("extractData Error: incompatible FunctionSpace");
2488     }
2489     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2490     if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2491     throw DataException("extractData Error: incompatible FunctionSpace");
2492     }
2493     }
2494     if (dataType==2) {
2495     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2496     if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2497     throw DataException("extractData Error: incompatible FunctionSpace");
2498     }
2499     }
2500     }
2501    
2502     //
2503     // Construct a DataVector to hold underlying data values
2504     DataVector dataVec(dataLength);
2505    
2506     //
2507     // Load this DataVector with the appropriate values
2508 jgs 123 int noValues;
2509     archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2510     cout << "\tnoValues: " << noValues << endl;
2511 jgs 119 switch (dataType) {
2512     case 0:
2513     // DataEmpty
2514 jgs 123 if (noValues != 0) {
2515     throw DataException("extractData Error: problem reading data from archive file");
2516     }
2517 jgs 119 break;
2518     case 1:
2519     // DataConstant
2520 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2521     throw DataException("extractData Error: problem reading data from archive file");
2522     }
2523 jgs 119 break;
2524     case 2:
2525     // DataTagged
2526 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2527     throw DataException("extractData Error: problem reading data from archive file");
2528     }
2529 jgs 119 break;
2530     case 3:
2531     // DataExpanded
2532 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2533     throw DataException("extractData Error: problem reading data from archive file");
2534     }
2535 jgs 119 break;
2536