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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26