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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 151 - (hide annotations)
Thu Sep 22 01:55:00 2005 UTC (14 years ago) by jgs
Original Path: trunk/esys2/escript/src/Data/Data.cpp
File size: 52056 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-22

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26