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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 782 - (hide annotations)
Tue Jul 18 00:47:47 2006 UTC (13 years, 8 months ago) by bcumming
File size: 60681 byte(s)
Large number of changes to Finley for meshing in MPI.

- optimisation and neatening up of rectcanglular mesh generation code
- first and second order 1D, 2D and 3D rectangular meshes are now
  available in finley and escript using MPI.
- reduced meshes now generated in MPI, and interpolation to and from 
  reduced data types now supported.  

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26