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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 149 - (hide annotations)
Thu Sep 1 03:31:39 2005 UTC (14 years, 1 month ago) by jgs
Original Path: trunk/esys2/escript/src/Data/Data.cpp
File size: 49419 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-01

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26