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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 286 - (hide annotations)
Thu Dec 1 06:40:20 2005 UTC (13 years, 11 months ago) by gross
Original Path: trunk/escript/src/Data/Data.cpp
File size: 52097 byte(s)
length method is removed as it is too slow. use length in util.py instead
1 jgs 94 // $Id$
2     /*=============================================================================
3    
4     ******************************************************************************
5     * *
6     * COPYRIGHT ACcESS 2004 - All Rights Reserved *
7     * *
8     * This software is the property of ACcESS. No part of this code *
9     * may be copied in any form or by any means without the expressed written *
10     * consent of ACcESS. Copying, use or modification of this software *
11     * by any unauthorised person is illegal unless that *
12     * person has a software license agreement with ACcESS. *
13     * *
14     ******************************************************************************
15    
16     ******************************************************************************/
17    
18     #include "escript/Data/Data.h"
19    
20     #include <iostream>
21 jgs 119 #include <fstream>
22 jgs 94 #include <algorithm>
23     #include <vector>
24     #include <exception>
25     #include <functional>
26     #include <math.h>
27    
28 jgs 153 #include <boost/python/dict.hpp>
29 jgs 94 #include <boost/python/str.hpp>
30     #include <boost/python/extract.hpp>
31     #include <boost/python/long.hpp>
32    
33     #include "escript/Data/DataException.h"
34     #include "escript/Data/DataExpanded.h"
35     #include "escript/Data/DataConstant.h"
36     #include "escript/Data/DataTagged.h"
37     #include "escript/Data/DataEmpty.h"
38     #include "escript/Data/DataArray.h"
39 jgs 123 #include "escript/Data/DataProf.h"
40 jgs 94 #include "escript/Data/FunctionSpaceFactory.h"
41     #include "escript/Data/AbstractContinuousDomain.h"
42 jgs 102 #include "escript/Data/UnaryFuncs.h"
43 jgs 94
44     using namespace std;
45     using namespace boost::python;
46     using namespace boost;
47     using namespace escript;
48    
49 jgs 151 #if defined DOPROF
50 jgs 123 //
51     // global table of profiling data for all Data objects
52     DataProf dataProfTable;
53 jgs 151 #endif
54 jgs 123
55 jgs 94 Data::Data()
56     {
57     //
58     // Default data is type DataEmpty
59     DataAbstract* temp=new DataEmpty();
60 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
61     m_data=temp_data;
62 jgs 151 #if defined DOPROF
63 jgs 123 // create entry in global profiling table for this object
64     profData = dataProfTable.newData();
65 jgs 151 #endif
66 jgs 94 }
67    
68     Data::Data(double value,
69     const tuple& shape,
70     const FunctionSpace& what,
71     bool expanded)
72     {
73     DataArrayView::ShapeType dataPointShape;
74     for (int i = 0; i < shape.attr("__len__")(); ++i) {
75     dataPointShape.push_back(extract<const int>(shape[i]));
76     }
77     DataArray temp(dataPointShape,value);
78     initialise(temp.getView(),what,expanded);
79 jgs 151 #if defined DOPROF
80 jgs 123 // create entry in global profiling table for this object
81     profData = dataProfTable.newData();
82 jgs 151 #endif
83 jgs 94 }
84    
85     Data::Data(double value,
86     const DataArrayView::ShapeType& dataPointShape,
87     const FunctionSpace& what,
88     bool expanded)
89     {
90     DataArray temp(dataPointShape,value);
91     pair<int,int> dataShape=what.getDataShape();
92     initialise(temp.getView(),what,expanded);
93 jgs 151 #if defined DOPROF
94 jgs 123 // create entry in global profiling table for this object
95     profData = dataProfTable.newData();
96 jgs 151 #endif
97 jgs 94 }
98    
99 jgs 102 Data::Data(const Data& inData)
100 jgs 94 {
101 jgs 102 m_data=inData.m_data;
102 jgs 151 #if defined DOPROF
103 jgs 123 // create entry in global profiling table for this object
104     profData = dataProfTable.newData();
105 jgs 151 #endif
106 jgs 94 }
107    
108     Data::Data(const Data& inData,
109     const DataArrayView::RegionType& region)
110     {
111     //
112 jgs 102 // Create Data which is a slice of another Data
113     DataAbstract* tmp = inData.m_data->getSlice(region);
114     shared_ptr<DataAbstract> temp_data(tmp);
115     m_data=temp_data;
116 jgs 151 #if defined DOPROF
117 jgs 123 // create entry in global profiling table for this object
118     profData = dataProfTable.newData();
119 jgs 151 #endif
120 jgs 94 }
121    
122     Data::Data(const Data& inData,
123     const FunctionSpace& functionspace)
124     {
125     if (inData.getFunctionSpace()==functionspace) {
126     m_data=inData.m_data;
127     } else {
128     Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
129 jgs 123 // Note: Must use a reference or pointer to a derived object
130 jgs 94 // in order to get polymorphic behaviour. Shouldn't really
131     // be able to create an instance of AbstractDomain but that was done
132 jgs 123 // as a boost:python work around which may no longer be required.
133 jgs 94 const AbstractDomain& inDataDomain=inData.getDomain();
134     if (inDataDomain==functionspace.getDomain()) {
135     inDataDomain.interpolateOnDomain(tmp,inData);
136     } else {
137     inDataDomain.interpolateACross(tmp,inData);
138     }
139     m_data=tmp.m_data;
140     }
141 jgs 151 #if defined DOPROF
142 jgs 123 // create entry in global profiling table for this object
143     profData = dataProfTable.newData();
144 jgs 151 #endif
145 jgs 94 }
146    
147     Data::Data(const DataTagged::TagListType& tagKeys,
148     const DataTagged::ValueListType & values,
149     const DataArrayView& defaultValue,
150     const FunctionSpace& what,
151     bool expanded)
152     {
153     DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
154 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
155     m_data=temp_data;
156 jgs 94 if (expanded) {
157     expand();
158     }
159 jgs 151 #if defined DOPROF
160 jgs 123 // create entry in global profiling table for this object
161     profData = dataProfTable.newData();
162 jgs 151 #endif
163 jgs 94 }
164    
165     Data::Data(const numeric::array& value,
166     const FunctionSpace& what,
167     bool expanded)
168     {
169     initialise(value,what,expanded);
170 jgs 151 #if defined DOPROF
171 jgs 123 // create entry in global profiling table for this object
172     profData = dataProfTable.newData();
173 jgs 151 #endif
174 jgs 94 }
175    
176     Data::Data(const DataArrayView& value,
177     const FunctionSpace& what,
178     bool expanded)
179     {
180     initialise(value,what,expanded);
181 jgs 151 #if defined DOPROF
182 jgs 123 // create entry in global profiling table for this object
183     profData = dataProfTable.newData();
184 jgs 151 #endif
185 jgs 94 }
186    
187     Data::Data(const object& value,
188     const FunctionSpace& what,
189     bool expanded)
190     {
191     numeric::array asNumArray(value);
192     initialise(asNumArray,what,expanded);
193 jgs 151 #if defined DOPROF
194 jgs 123 // create entry in global profiling table for this object
195     profData = dataProfTable.newData();
196 jgs 151 #endif
197 jgs 94 }
198    
199     Data::Data(const object& value,
200     const Data& other)
201     {
202     //
203     // Create DataConstant using the given value and all other parameters
204     // copied from other. If value is a rank 0 object this Data
205     // will assume the point data shape of other.
206     DataArray temp(value);
207     if (temp.getView().getRank()==0) {
208     //
209     // Create a DataArray with the scalar value for all elements
210     DataArray temp2(other.getPointDataView().getShape(),temp.getView()());
211     initialise(temp2.getView(),other.getFunctionSpace(),false);
212     } else {
213     //
214     // Create a DataConstant with the same sample shape as other
215     initialise(temp.getView(),other.getFunctionSpace(),false);
216     }
217 jgs 151 #if defined DOPROF
218 jgs 123 // create entry in global profiling table for this object
219     profData = dataProfTable.newData();
220 jgs 151 #endif
221 jgs 94 }
222    
223 jgs 151 Data::~Data()
224     {
225    
226     }
227    
228 jgs 94 escriptDataC
229     Data::getDataC()
230     {
231     escriptDataC temp;
232     temp.m_dataPtr=(void*)this;
233     return temp;
234     }
235    
236     escriptDataC
237     Data::getDataC() const
238     {
239     escriptDataC temp;
240     temp.m_dataPtr=(void*)this;
241     return temp;
242     }
243    
244 jgs 121 const boost::python::tuple
245 jgs 94 Data::getShapeTuple() const
246     {
247     const DataArrayView::ShapeType& shape=getDataPointShape();
248     switch(getDataPointRank()) {
249     case 0:
250     return make_tuple();
251     case 1:
252     return make_tuple(long_(shape[0]));
253     case 2:
254     return make_tuple(long_(shape[0]),long_(shape[1]));
255     case 3:
256     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]));
257     case 4:
258     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]),long_(shape[3]));
259     default:
260     throw DataException("Error - illegal Data rank.");
261     }
262     }
263    
264     void
265     Data::copy(const Data& other)
266     {
267     //
268     // Perform a deep copy
269     {
270     DataExpanded* temp=dynamic_cast<DataExpanded*>(other.m_data.get());
271     if (temp!=0) {
272     //
273     // Construct a DataExpanded copy
274     DataAbstract* newData=new DataExpanded(*temp);
275 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
276     m_data=temp_data;
277 jgs 94 return;
278     }
279     }
280     {
281     DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());
282     if (temp!=0) {
283     //
284 jgs 102 // Construct a DataTagged copy
285 jgs 94 DataAbstract* newData=new DataTagged(*temp);
286 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
287     m_data=temp_data;
288 jgs 94 return;
289     }
290     }
291     {
292     DataConstant* temp=dynamic_cast<DataConstant*>(other.m_data.get());
293     if (temp!=0) {
294     //
295     // Construct a DataConstant copy
296     DataAbstract* newData=new DataConstant(*temp);
297 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
298     m_data=temp_data;
299 jgs 94 return;
300     }
301     }
302 jgs 102 {
303     DataEmpty* temp=dynamic_cast<DataEmpty*>(other.m_data.get());
304     if (temp!=0) {
305     //
306     // Construct a DataEmpty copy
307     DataAbstract* newData=new DataEmpty();
308     shared_ptr<DataAbstract> temp_data(newData);
309     m_data=temp_data;
310     return;
311     }
312     }
313 jgs 94 throw DataException("Error - Copy not implemented for this Data type.");
314     }
315    
316     void
317     Data::copyWithMask(const Data& other,
318     const Data& mask)
319     {
320     Data mask1;
321     Data mask2;
322    
323     mask1 = mask.wherePositive();
324     mask2.copy(mask1);
325    
326     mask1 *= other;
327     mask2 *= *this;
328     mask2 = *this - mask2;
329    
330     *this = mask1 + mask2;
331     }
332    
333     bool
334     Data::isExpanded() const
335     {
336     DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
337     return (temp!=0);
338     }
339    
340     bool
341     Data::isTagged() const
342     {
343     DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
344     return (temp!=0);
345     }
346    
347     bool
348     Data::isEmpty() const
349     {
350     DataEmpty* temp=dynamic_cast<DataEmpty*>(m_data.get());
351     return (temp!=0);
352     }
353    
354     bool
355     Data::isConstant() const
356     {
357     DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
358     return (temp!=0);
359     }
360    
361     void
362     Data::expand()
363     {
364     if (isConstant()) {
365     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
366     DataAbstract* temp=new DataExpanded(*tempDataConst);
367 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
368     m_data=temp_data;
369 jgs 94 } else if (isTagged()) {
370     DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
371     DataAbstract* temp=new DataExpanded(*tempDataTag);
372 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
373     m_data=temp_data;
374 jgs 94 } else if (isExpanded()) {
375     //
376     // do nothing
377     } else if (isEmpty()) {
378     throw DataException("Error - Expansion of DataEmpty not possible.");
379     } else {
380     throw DataException("Error - Expansion not implemented for this Data type.");
381     }
382     }
383    
384     void
385     Data::tag()
386     {
387     if (isConstant()) {
388     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
389     DataAbstract* temp=new DataTagged(*tempDataConst);
390 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
391     m_data=temp_data;
392 jgs 94 } else if (isTagged()) {
393     // do nothing
394     } else if (isExpanded()) {
395     throw DataException("Error - Creating tag data from DataExpanded not possible.");
396     } else if (isEmpty()) {
397     throw DataException("Error - Creating tag data from DataEmpty not possible.");
398     } else {
399     throw DataException("Error - Tagging not implemented for this Data type.");
400     }
401     }
402    
403 jgs 102 void
404     Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)
405     {
406     m_data->reshapeDataPoint(shape);
407     }
408    
409 jgs 94 Data
410     Data::wherePositive() const
411     {
412 jgs 151 #if defined DOPROF
413 jgs 123 profData->where++;
414 jgs 151 #endif
415 jgs 94 return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
416     }
417    
418     Data
419 jgs 102 Data::whereNegative() const
420     {
421 jgs 151 #if defined DOPROF
422 jgs 123 profData->where++;
423 jgs 151 #endif
424 jgs 102 return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
425     }
426    
427     Data
428 jgs 94 Data::whereNonNegative() const
429     {
430 jgs 151 #if defined DOPROF
431 jgs 123 profData->where++;
432 jgs 151 #endif
433 jgs 94 return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
434     }
435    
436     Data
437 jgs 102 Data::whereNonPositive() const
438 jgs 94 {
439 jgs 151 #if defined DOPROF
440 jgs 123 profData->where++;
441 jgs 151 #endif
442 jgs 102 return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
443 jgs 94 }
444    
445     Data
446     Data::whereZero() const
447     {
448 jgs 151 #if defined DOPROF
449 jgs 123 profData->where++;
450 jgs 151 #endif
451 jgs 94 return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));
452     }
453    
454     Data
455 jgs 102 Data::whereNonZero() const
456     {
457 jgs 151 #if defined DOPROF
458 jgs 123 profData->where++;
459 jgs 151 #endif
460 jgs 102 return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));
461     }
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     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[i,j] = integrals[index];
888     }
889     }
890     }
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     bp_array[i,j,k] = integrals[index];
898     }
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     bp_array[i,j,k,l] = integrals[index];
910     }
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