/[escript]/branches/schroedinger/escript/src/Data.cpp
ViewVC logotype

Annotation of /branches/schroedinger/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 480 - (hide annotations)
Wed Feb 1 05:15:12 2006 UTC (13 years, 3 months ago) by jgs
Original Path: trunk/escript/src/Data.cpp
File size: 51939 byte(s)
rationalise #includes and forward declarations

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26