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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 153 - (hide annotations)
Tue Oct 25 01:51:20 2005 UTC (13 years, 8 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/Data.cpp
File size: 52231 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-10-25

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 jgs 94 Data::log() const
1031     {
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     Data::ln() const
1040     {
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::length() const
1178     {
1179 jgs 151 #if defined DOPROF
1180 jgs 123 profData->reduction2++;
1181 jgs 151 #endif
1182 jgs 147 Length len_func;
1183     return dp_algorithm(len_func,0);
1184 jgs 123 }
1185    
1186     Data
1187     Data::trace() const
1188     {
1189 jgs 151 #if defined DOPROF
1190 jgs 123 profData->reduction2++;
1191 jgs 151 #endif
1192 jgs 147 Trace trace_func;
1193     return dp_algorithm(trace_func,0);
1194 jgs 123 }
1195    
1196     Data
1197     Data::transpose(int axis) const
1198     {
1199 jgs 151 #if defined DOPROF
1200 jgs 123 profData->reduction2++;
1201 jgs 151 #endif
1202 jgs 123 // not implemented
1203     throw DataException("Error - Data::transpose not implemented yet.");
1204     return Data();
1205     }
1206    
1207 jgs 121 const boost::python::tuple
1208     Data::mindp() const
1209     {
1210 jgs 148 // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1211     // abort (for unknown reasons) if there are openmp directives with it in the
1212     // surrounding function
1213    
1214     int SampleNo;
1215     int DataPointNo;
1216    
1217     calc_mindp(SampleNo,DataPointNo);
1218    
1219     return make_tuple(SampleNo,DataPointNo);
1220     }
1221    
1222     void
1223     Data::calc_mindp(int& SampleNo,
1224     int& DataPointNo) const
1225     {
1226     int i,j;
1227     int lowi=0,lowj=0;
1228     double min=numeric_limits<double>::max();
1229    
1230 jgs 121 Data temp=minval();
1231    
1232     int numSamples=temp.getNumSamples();
1233     int numDPPSample=temp.getNumDataPointsPerSample();
1234    
1235 jgs 148 double next,local_min;
1236     int local_lowi,local_lowj;
1237 jgs 121
1238 jgs 148 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1239     {
1240     local_min=min;
1241     #pragma omp for private(i,j) schedule(static)
1242     for (i=0; i<numSamples; i++) {
1243     for (j=0; j<numDPPSample; j++) {
1244     next=temp.getDataPoint(i,j)();
1245     if (next<local_min) {
1246     local_min=next;
1247     local_lowi=i;
1248     local_lowj=j;
1249     }
1250 jgs 121 }
1251     }
1252 jgs 148 #pragma omp critical
1253     if (local_min<min) {
1254     min=local_min;
1255     lowi=local_lowi;
1256     lowj=local_lowj;
1257     }
1258 jgs 121 }
1259    
1260 jgs 148 SampleNo = lowi;
1261     DataPointNo = lowj;
1262 jgs 121 }
1263    
1264 jgs 104 void
1265     Data::saveDX(std::string fileName) const
1266     {
1267 jgs 153 boost::python::dict args;
1268     args["data"]=boost::python::object(this);
1269     getDomain().saveDX(fileName,args);
1270 jgs 104 return;
1271     }
1272    
1273 jgs 110 void
1274     Data::saveVTK(std::string fileName) const
1275     {
1276 jgs 153 boost::python::dict args;
1277     args["data"]=boost::python::object(this);
1278     getDomain().saveVTK(fileName,args);
1279 jgs 110 return;
1280     }
1281    
1282 jgs 102 Data&
1283     Data::operator+=(const Data& right)
1284     {
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 boost::python::object& right)
1294 jgs 94 {
1295 jgs 151 #if defined DOPROF
1296 jgs 123 profData->binary++;
1297 jgs 151 #endif
1298 jgs 94 binaryOp(right,plus<double>());
1299     return (*this);
1300     }
1301    
1302 jgs 102 Data&
1303     Data::operator-=(const Data& 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 boost::python::object& right)
1314 jgs 94 {
1315 jgs 151 #if defined DOPROF
1316 jgs 123 profData->binary++;
1317 jgs 151 #endif
1318 jgs 94 binaryOp(right,minus<double>());
1319     return (*this);
1320     }
1321    
1322 jgs 102 Data&
1323     Data::operator*=(const Data& 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 boost::python::object& right)
1334 jgs 94 {
1335 jgs 151 #if defined DOPROF
1336 jgs 123 profData->binary++;
1337 jgs 151 #endif
1338 jgs 94 binaryOp(right,multiplies<double>());
1339     return (*this);
1340     }
1341    
1342 jgs 102 Data&
1343     Data::operator/=(const Data& 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::operator/=(const boost::python::object& right)
1354 jgs 94 {
1355 jgs 151 #if defined DOPROF
1356 jgs 123 profData->binary++;
1357 jgs 151 #endif
1358 jgs 94 binaryOp(right,divides<double>());
1359     return (*this);
1360     }
1361    
1362 jgs 102 Data
1363     Data::powO(const boost::python::object& right) const
1364 jgs 94 {
1365 jgs 151 #if defined DOPROF
1366 jgs 123 profData->binary++;
1367 jgs 151 #endif
1368 jgs 94 Data result;
1369     result.copy(*this);
1370     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1371     return result;
1372     }
1373    
1374 jgs 102 Data
1375     Data::powD(const Data& right) const
1376 jgs 94 {
1377 jgs 151 #if defined DOPROF
1378 jgs 123 profData->binary++;
1379 jgs 151 #endif
1380 jgs 94 Data result;
1381     result.copy(*this);
1382     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1383     return result;
1384     }
1385    
1386     //
1387 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1388 jgs 102 Data
1389     escript::operator+(const Data& left, const Data& right)
1390 jgs 94 {
1391     Data result;
1392     //
1393     // perform a deep copy
1394     result.copy(left);
1395     result+=right;
1396     return result;
1397     }
1398    
1399     //
1400 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1401 jgs 102 Data
1402     escript::operator-(const Data& left, const Data& right)
1403 jgs 94 {
1404     Data result;
1405     //
1406     // perform a deep copy
1407     result.copy(left);
1408     result-=right;
1409     return result;
1410     }
1411    
1412     //
1413 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1414 jgs 102 Data
1415     escript::operator*(const Data& left, const Data& right)
1416 jgs 94 {
1417     Data result;
1418     //
1419     // perform a deep copy
1420     result.copy(left);
1421     result*=right;
1422     return result;
1423     }
1424    
1425     //
1426 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1427 jgs 102 Data
1428     escript::operator/(const Data& left, const Data& right)
1429 jgs 94 {
1430     Data result;
1431     //
1432     // perform a deep copy
1433     result.copy(left);
1434     result/=right;
1435     return result;
1436     }
1437    
1438     //
1439 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1440 jgs 102 Data
1441     escript::operator+(const Data& left, const boost::python::object& right)
1442 jgs 94 {
1443     //
1444     // Convert to DataArray format if possible
1445     DataArray temp(right);
1446     Data result;
1447     //
1448     // perform a deep copy
1449     result.copy(left);
1450     result+=right;
1451     return result;
1452     }
1453    
1454     //
1455 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1456 jgs 102 Data
1457     escript::operator-(const Data& left, const boost::python::object& right)
1458 jgs 94 {
1459     //
1460     // Convert to DataArray format if possible
1461     DataArray temp(right);
1462     Data result;
1463     //
1464     // perform a deep copy
1465     result.copy(left);
1466     result-=right;
1467     return result;
1468     }
1469    
1470     //
1471 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1472 jgs 102 Data
1473     escript::operator*(const Data& left, const boost::python::object& right)
1474 jgs 94 {
1475     //
1476     // Convert to DataArray format if possible
1477     DataArray temp(right);
1478     Data result;
1479     //
1480     // perform a deep copy
1481     result.copy(left);
1482     result*=right;
1483     return result;
1484     }
1485    
1486     //
1487 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1488 jgs 102 Data
1489     escript::operator/(const Data& left, const boost::python::object& right)
1490 jgs 94 {
1491     //
1492     // Convert to DataArray format if possible
1493     DataArray temp(right);
1494     Data result;
1495     //
1496     // perform a deep copy
1497     result.copy(left);
1498     result/=right;
1499     return result;
1500     }
1501    
1502     //
1503 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1504 jgs 102 Data
1505     escript::operator+(const boost::python::object& left, const Data& right)
1506 jgs 94 {
1507     //
1508     // Construct the result using the given value and the other parameters
1509     // from right
1510     Data result(left,right);
1511     result+=right;
1512     return result;
1513     }
1514    
1515     //
1516 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1517 jgs 102 Data
1518     escript::operator-(const boost::python::object& left, const Data& right)
1519 jgs 94 {
1520     //
1521     // Construct the result using the given value and the other parameters
1522     // from right
1523     Data result(left,right);
1524     result-=right;
1525     return result;
1526     }
1527    
1528     //
1529 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1530 jgs 102 Data
1531     escript::operator*(const boost::python::object& left, const Data& right)
1532 jgs 94 {
1533     //
1534     // Construct the result using the given value and the other parameters
1535     // from right
1536     Data result(left,right);
1537     result*=right;
1538     return result;
1539     }
1540    
1541     //
1542 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1543 jgs 102 Data
1544     escript::operator/(const boost::python::object& left, const Data& right)
1545 jgs 94 {
1546     //
1547     // Construct the result using the given value and the other parameters
1548     // from right
1549     Data result(left,right);
1550     result/=right;
1551     return result;
1552     }
1553    
1554     //
1555 jgs 102 //bool escript::operator==(const Data& left, const Data& right)
1556     //{
1557     // /*
1558     // NB: this operator does very little at this point, and isn't to
1559     // be relied on. Requires further implementation.
1560     // */
1561     //
1562     // bool ret;
1563     //
1564     // if (left.isEmpty()) {
1565     // if(!right.isEmpty()) {
1566     // ret = false;
1567     // } else {
1568     // ret = true;
1569     // }
1570     // }
1571     //
1572     // if (left.isConstant()) {
1573     // if(!right.isConstant()) {
1574     // ret = false;
1575     // } else {
1576     // ret = true;
1577     // }
1578     // }
1579     //
1580     // if (left.isTagged()) {
1581     // if(!right.isTagged()) {
1582     // ret = false;
1583     // } else {
1584     // ret = true;
1585     // }
1586     // }
1587     //
1588     // if (left.isExpanded()) {
1589     // if(!right.isExpanded()) {
1590     // ret = false;
1591     // } else {
1592     // ret = true;
1593     // }
1594     // }
1595     //
1596     // return ret;
1597     //}
1598    
1599     Data
1600     Data::getItem(const boost::python::object& key) const
1601 jgs 94 {
1602 jgs 102 const DataArrayView& view=getPointDataView();
1603 jgs 94
1604 jgs 102 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1605 jgs 94
1606 jgs 102 if (slice_region.size()!=view.getRank()) {
1607     throw DataException("Error - slice size does not match Data rank.");
1608 jgs 94 }
1609    
1610 jgs 102 return getSlice(slice_region);
1611 jgs 94 }
1612    
1613     Data
1614 jgs 102 Data::getSlice(const DataArrayView::RegionType& region) const
1615 jgs 94 {
1616 jgs 151 #if defined DOPROF
1617 jgs 123 profData->slicing++;
1618 jgs 151 #endif
1619 jgs 102 return Data(*this,region);
1620 jgs 94 }
1621    
1622     void
1623 jgs 102 Data::setItemO(const boost::python::object& key,
1624     const boost::python::object& value)
1625 jgs 94 {
1626 jgs 102 Data tempData(value,getFunctionSpace());
1627     setItemD(key,tempData);
1628     }
1629    
1630     void
1631     Data::setItemD(const boost::python::object& key,
1632     const Data& value)
1633     {
1634 jgs 94 const DataArrayView& view=getPointDataView();
1635 jgs 123
1636 jgs 94 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1637     if (slice_region.size()!=view.getRank()) {
1638     throw DataException("Error - slice size does not match Data rank.");
1639     }
1640 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
1641     setSlice(Data(value,getFunctionSpace()),slice_region);
1642     } else {
1643     setSlice(value,slice_region);
1644     }
1645 jgs 94 }
1646    
1647     void
1648 jgs 102 Data::setSlice(const Data& value,
1649     const DataArrayView::RegionType& region)
1650 jgs 94 {
1651 jgs 151 #if defined DOPROF
1652 jgs 123 profData->slicing++;
1653 jgs 151 #endif
1654 jgs 102 Data tempValue(value);
1655     typeMatchLeft(tempValue);
1656     typeMatchRight(tempValue);
1657     m_data->setSlice(tempValue.m_data.get(),region);
1658     }
1659    
1660     void
1661     Data::typeMatchLeft(Data& right) const
1662     {
1663     if (isExpanded()){
1664     right.expand();
1665     } else if (isTagged()) {
1666     if (right.isConstant()) {
1667     right.tag();
1668     }
1669     }
1670     }
1671    
1672     void
1673     Data::typeMatchRight(const Data& right)
1674     {
1675 jgs 94 if (isTagged()) {
1676     if (right.isExpanded()) {
1677     expand();
1678     }
1679     } else if (isConstant()) {
1680     if (right.isExpanded()) {
1681     expand();
1682     } else if (right.isTagged()) {
1683     tag();
1684     }
1685     }
1686     }
1687    
1688     void
1689     Data::setTaggedValue(int tagKey,
1690     const boost::python::object& value)
1691     {
1692     //
1693     // Ensure underlying data object is of type DataTagged
1694     tag();
1695    
1696     if (!isTagged()) {
1697     throw DataException("Error - DataTagged conversion failed!!");
1698     }
1699    
1700     //
1701     // Construct DataArray from boost::python::object input value
1702     DataArray valueDataArray(value);
1703    
1704     //
1705     // Call DataAbstract::setTaggedValue
1706     m_data->setTaggedValue(tagKey,valueDataArray.getView());
1707     }
1708    
1709 jgs 110 void
1710 jgs 121 Data::setTaggedValueFromCPP(int tagKey,
1711     const DataArrayView& value)
1712     {
1713     //
1714     // Ensure underlying data object is of type DataTagged
1715     tag();
1716    
1717     if (!isTagged()) {
1718     throw DataException("Error - DataTagged conversion failed!!");
1719     }
1720    
1721     //
1722     // Call DataAbstract::setTaggedValue
1723     m_data->setTaggedValue(tagKey,value);
1724     }
1725    
1726 jgs 149 int
1727     Data::getTagNumber(int dpno)
1728     {
1729     return m_data->getTagNumber(dpno);
1730     }
1731    
1732 jgs 121 void
1733 jgs 110 Data::setRefValue(int ref,
1734     const boost::python::numeric::array& value)
1735     {
1736     //
1737     // Construct DataArray from boost::python::object input value
1738     DataArray valueDataArray(value);
1739    
1740     //
1741     // Call DataAbstract::setRefValue
1742     m_data->setRefValue(ref,valueDataArray);
1743     }
1744    
1745     void
1746     Data::getRefValue(int ref,
1747     boost::python::numeric::array& value)
1748     {
1749     //
1750     // Construct DataArray for boost::python::object return value
1751     DataArray valueDataArray(value);
1752    
1753     //
1754     // Load DataArray with values from data-points specified by ref
1755     m_data->getRefValue(ref,valueDataArray);
1756    
1757     //
1758     // Load values from valueDataArray into return numarray
1759    
1760     // extract the shape of the numarray
1761     int rank = value.getrank();
1762     DataArrayView::ShapeType shape;
1763     for (int i=0; i < rank; i++) {
1764     shape.push_back(extract<int>(value.getshape()[i]));
1765     }
1766    
1767     // and load the numarray with the data from the DataArray
1768     DataArrayView valueView = valueDataArray.getView();
1769    
1770     if (rank==0) {
1771 jgs 126 boost::python::numeric::array temp_numArray(valueView());
1772     value = temp_numArray;
1773 jgs 110 }
1774     if (rank==1) {
1775     for (int i=0; i < shape[0]; i++) {
1776     value[i] = valueView(i);
1777     }
1778     }
1779     if (rank==2) {
1780 jgs 126 for (int i=0; i < shape[0]; i++) {
1781     for (int j=0; j < shape[1]; j++) {
1782     value[i][j] = valueView(i,j);
1783     }
1784     }
1785 jgs 110 }
1786     if (rank==3) {
1787 jgs 126 for (int i=0; i < shape[0]; i++) {
1788     for (int j=0; j < shape[1]; j++) {
1789     for (int k=0; k < shape[2]; k++) {
1790     value[i][j][k] = valueView(i,j,k);
1791     }
1792     }
1793     }
1794 jgs 110 }
1795     if (rank==4) {
1796 jgs 126 for (int i=0; i < shape[0]; i++) {
1797     for (int j=0; j < shape[1]; j++) {
1798     for (int k=0; k < shape[2]; k++) {
1799     for (int l=0; l < shape[3]; l++) {
1800     value[i][j][k][l] = valueView(i,j,k,l);
1801     }
1802     }
1803     }
1804     }
1805 jgs 110 }
1806    
1807     }
1808    
1809 jgs 94 void
1810 jgs 119 Data::archiveData(const std::string fileName)
1811     {
1812     cout << "Archiving Data object to: " << fileName << endl;
1813    
1814     //
1815     // Determine type of this Data object
1816     int dataType = -1;
1817    
1818     if (isEmpty()) {
1819     dataType = 0;
1820     cout << "\tdataType: DataEmpty" << endl;
1821     }
1822     if (isConstant()) {
1823     dataType = 1;
1824     cout << "\tdataType: DataConstant" << endl;
1825     }
1826     if (isTagged()) {
1827     dataType = 2;
1828     cout << "\tdataType: DataTagged" << endl;
1829     }
1830     if (isExpanded()) {
1831     dataType = 3;
1832     cout << "\tdataType: DataExpanded" << endl;
1833     }
1834 jgs 123
1835 jgs 119 if (dataType == -1) {
1836     throw DataException("archiveData Error: undefined dataType");
1837     }
1838    
1839     //
1840     // Collect data items common to all Data types
1841     int noSamples = getNumSamples();
1842     int noDPPSample = getNumDataPointsPerSample();
1843     int functionSpaceType = getFunctionSpace().getTypeCode();
1844     int dataPointRank = getDataPointRank();
1845     int dataPointSize = getDataPointSize();
1846     int dataLength = getLength();
1847     DataArrayView::ShapeType dataPointShape = getDataPointShape();
1848     int referenceNumbers[noSamples];
1849     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1850     referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1851     }
1852     int tagNumbers[noSamples];
1853     if (isTagged()) {
1854     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1855     tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1856     }
1857     }
1858    
1859     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1860     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1861     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1862    
1863     //
1864     // Flatten Shape to an array of integers suitable for writing to file
1865     int flatShape[4] = {0,0,0,0};
1866     cout << "\tshape: < ";
1867     for (int dim=0; dim<dataPointRank; dim++) {
1868     flatShape[dim] = dataPointShape[dim];
1869     cout << dataPointShape[dim] << " ";
1870     }
1871     cout << ">" << endl;
1872    
1873     //
1874 jgs 123 // Open archive file
1875 jgs 119 ofstream archiveFile;
1876     archiveFile.open(fileName.data(), ios::out);
1877    
1878     if (!archiveFile.good()) {
1879     throw DataException("archiveData Error: problem opening archive file");
1880     }
1881    
1882 jgs 123 //
1883     // Write common data items to archive file
1884 jgs 119 archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1885     archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1886     archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1887     archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1888     archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1889     archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1890     archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1891     for (int dim = 0; dim < 4; dim++) {
1892     archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1893     }
1894     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1895     archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1896     }
1897     if (isTagged()) {
1898     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1899     archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1900     }
1901     }
1902    
1903     if (!archiveFile.good()) {
1904     throw DataException("archiveData Error: problem writing to archive file");
1905     }
1906    
1907     //
1908 jgs 123 // Archive underlying data values for each Data type
1909     int noValues;
1910 jgs 119 switch (dataType) {
1911     case 0:
1912     // DataEmpty
1913 jgs 123 noValues = 0;
1914     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1915     cout << "\tnoValues: " << noValues << endl;
1916 jgs 119 break;
1917     case 1:
1918     // DataConstant
1919 jgs 123 noValues = m_data->getLength();
1920     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1921     cout << "\tnoValues: " << noValues << endl;
1922     if (m_data->archiveData(archiveFile,noValues)) {
1923     throw DataException("archiveData Error: problem writing data to archive file");
1924     }
1925 jgs 119 break;
1926     case 2:
1927     // DataTagged
1928 jgs 123 noValues = m_data->getLength();
1929     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1930     cout << "\tnoValues: " << noValues << endl;
1931     if (m_data->archiveData(archiveFile,noValues)) {
1932     throw DataException("archiveData Error: problem writing data to archive file");
1933     }
1934 jgs 119 break;
1935     case 3:
1936     // DataExpanded
1937 jgs 123 noValues = m_data->getLength();
1938     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1939     cout << "\tnoValues: " << noValues << endl;
1940     if (m_data->archiveData(archiveFile,noValues)) {
1941     throw DataException("archiveData Error: problem writing data to archive file");
1942     }
1943 jgs 119 break;
1944     }
1945    
1946 jgs 123 if (!archiveFile.good()) {
1947     throw DataException("archiveData Error: problem writing data to archive file");
1948     }
1949    
1950     //
1951     // Close archive file
1952     archiveFile.close();
1953    
1954     if (!archiveFile.good()) {
1955     throw DataException("archiveData Error: problem closing archive file");
1956     }
1957    
1958 jgs 119 }
1959    
1960     void
1961     Data::extractData(const std::string fileName,
1962     const FunctionSpace& fspace)
1963     {
1964     //
1965     // Can only extract Data to an object which is initially DataEmpty
1966     if (!isEmpty()) {
1967     throw DataException("extractData Error: can only extract to DataEmpty object");
1968     }
1969    
1970     cout << "Extracting Data object from: " << fileName << endl;
1971    
1972     int dataType;
1973     int noSamples;
1974     int noDPPSample;
1975     int functionSpaceType;
1976     int dataPointRank;
1977     int dataPointSize;
1978     int dataLength;
1979     DataArrayView::ShapeType dataPointShape;
1980     int flatShape[4];
1981    
1982     //
1983 jgs 123 // Open the archive file
1984 jgs 119 ifstream archiveFile;
1985     archiveFile.open(fileName.data(), ios::in);
1986    
1987     if (!archiveFile.good()) {
1988     throw DataException("extractData Error: problem opening archive file");
1989     }
1990    
1991 jgs 123 //
1992     // Read common data items from archive file
1993 jgs 119 archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1994     archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1995     archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1996     archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1997     archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1998     archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1999     archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
2000     for (int dim = 0; dim < 4; dim++) {
2001     archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2002     if (flatShape[dim]>0) {
2003     dataPointShape.push_back(flatShape[dim]);
2004     }
2005     }
2006     int referenceNumbers[noSamples];
2007     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2008     archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2009     }
2010     int tagNumbers[noSamples];
2011     if (dataType==2) {
2012     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2013     archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2014     }
2015     }
2016    
2017     if (!archiveFile.good()) {
2018     throw DataException("extractData Error: problem reading from archive file");
2019     }
2020    
2021 jgs 123 //
2022     // Verify the values just read from the archive file
2023 jgs 119 switch (dataType) {
2024     case 0:
2025     cout << "\tdataType: DataEmpty" << endl;
2026     break;
2027     case 1:
2028     cout << "\tdataType: DataConstant" << endl;
2029     break;
2030     case 2:
2031     cout << "\tdataType: DataTagged" << endl;
2032     break;
2033     case 3:
2034     cout << "\tdataType: DataExpanded" << endl;
2035     break;
2036     default:
2037     throw DataException("extractData Error: undefined dataType read from archive file");
2038     break;
2039     }
2040    
2041     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2042     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2043     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2044     cout << "\tshape: < ";
2045     for (int dim = 0; dim < dataPointRank; dim++) {
2046     cout << dataPointShape[dim] << " ";
2047     }
2048     cout << ">" << endl;
2049    
2050     //
2051     // Verify that supplied FunctionSpace object is compatible with this Data object.
2052     if ( (fspace.getTypeCode()!=functionSpaceType) ||
2053     (fspace.getNumSamples()!=noSamples) ||
2054     (fspace.getNumDPPSample()!=noDPPSample)
2055     ) {
2056     throw DataException("extractData Error: incompatible FunctionSpace");
2057     }
2058     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2059     if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2060     throw DataException("extractData Error: incompatible FunctionSpace");
2061     }
2062     }
2063     if (dataType==2) {
2064     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2065     if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2066     throw DataException("extractData Error: incompatible FunctionSpace");
2067     }
2068     }
2069     }
2070    
2071     //
2072     // Construct a DataVector to hold underlying data values
2073     DataVector dataVec(dataLength);
2074    
2075     //
2076     // Load this DataVector with the appropriate values
2077 jgs 123 int noValues;
2078     archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2079     cout << "\tnoValues: " << noValues << endl;
2080 jgs 119 switch (dataType) {
2081     case 0:
2082     // DataEmpty
2083 jgs 123 if (noValues != 0) {
2084     throw DataException("extractData Error: problem reading data from archive file");
2085     }
2086 jgs 119 break;
2087     case 1:
2088     // DataConstant
2089 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2090     throw DataException("extractData Error: problem reading data from archive file");
2091     }
2092 jgs 119 break;
2093     case 2:
2094     // DataTagged
2095 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2096     throw DataException("extractData Error: problem reading data from archive file");
2097     }
2098 jgs 119 break;
2099     case 3:
2100     // DataExpanded
2101 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2102     throw DataException("extractData Error: problem reading data from archive file");
2103     }
2104 jgs 119 break;
2105     }
2106    
2107 jgs 123 if (!archiveFile.good()) {
2108     throw DataException("extractData Error: problem reading from archive file");
2109     }
2110    
2111 jgs 119 //
2112 jgs 123 // Close archive file
2113     archiveFile.close();
2114    
2115     if (!archiveFile.good()) {
2116     throw DataException("extractData Error: problem closing archive file");
2117     }
2118    
2119     //
2120 jgs 119 // Construct an appropriate Data object
2121     DataAbstract* tempData;
2122     switch (dataType) {
2123     case 0:
2124     // DataEmpty
2125     tempData=new DataEmpty();
2126     break;
2127     case 1:
2128     // DataConstant
2129     tempData=new DataConstant(fspace,dataPointShape,dataVec);
2130     break;
2131     case 2:
2132     // DataTagged
2133     tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2134     break;
2135     case 3:
2136     // DataExpanded
2137     tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2138     break;
2139     }
2140     shared_ptr<DataAbstract> temp_data(tempData);
2141     m_data=temp_data;
2142     }
2143    
2144 jgs 94 ostream& escript::operator<<(ostream& o, const Data& data)
2145     {
2146     o << data.toString();
2147     return o;
2148     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26