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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 474 - (hide annotations)
Mon Jan 30 04:23:44 2006 UTC (13 years, 10 months ago) by jgs
File size: 52045 byte(s)
restructure escript source tree
move src/Data/* -> src
remove inc
modify #includes and cpppath settings accordingly

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26