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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (13 years, 9 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/Data.cpp
File size: 50424 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

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 jgs 150 Data::asin() const
896     {
897     profData->unary++;
898     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
899     }
900    
901     Data
902     Data::acos() const
903     {
904     profData->unary++;
905     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
906     }
907    
908     Data
909     Data::atan() const
910     {
911     profData->unary++;
912     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
913     }
914    
915     Data
916     Data::sinh() const
917     {
918     profData->unary++;
919     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
920     }
921    
922     Data
923     Data::cosh() const
924     {
925     profData->unary++;
926     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
927     }
928    
929     Data
930     Data::tanh() const
931     {
932     profData->unary++;
933     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
934     }
935    
936     Data
937     Data::asinh() const
938     {
939     profData->unary++;
940     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
941     }
942    
943     Data
944     Data::acosh() const
945     {
946     profData->unary++;
947     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
948     }
949    
950     Data
951     Data::atanh() const
952     {
953     profData->unary++;
954     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
955     }
956    
957     Data
958 jgs 94 Data::log() const
959     {
960 jgs 123 profData->unary++;
961 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
962     }
963    
964     Data
965     Data::ln() const
966     {
967 jgs 123 profData->unary++;
968 jgs 94 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
969     }
970    
971 jgs 106 Data
972     Data::sign() const
973 jgs 94 {
974 jgs 123 profData->unary++;
975 jgs 106 return escript::unaryOp(*this,escript::fsign);
976 jgs 94 }
977    
978 jgs 106 Data
979     Data::abs() const
980 jgs 94 {
981 jgs 123 profData->unary++;
982 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
983 jgs 94 }
984    
985 jgs 106 Data
986     Data::neg() const
987 jgs 94 {
988 jgs 123 profData->unary++;
989 jgs 106 return escript::unaryOp(*this,negate<double>());
990 jgs 94 }
991    
992 jgs 102 Data
993 jgs 106 Data::pos() const
994 jgs 94 {
995 jgs 123 profData->unary++;
996 jgs 148 Data result;
997     // perform a deep copy
998     result.copy(*this);
999     return result;
1000 jgs 102 }
1001    
1002     Data
1003 jgs 106 Data::exp() const
1004 jgs 102 {
1005 jgs 123 profData->unary++;
1006 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1007 jgs 102 }
1008    
1009     Data
1010 jgs 106 Data::sqrt() const
1011 jgs 102 {
1012 jgs 123 profData->unary++;
1013 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1014 jgs 102 }
1015    
1016 jgs 106 double
1017     Data::Lsup() const
1018 jgs 102 {
1019 jgs 123 profData->reduction1++;
1020 jgs 106 //
1021     // set the initial absolute maximum value to zero
1022 jgs 147 AbsMax abs_max_func;
1023     return algorithm(abs_max_func,0);
1024 jgs 102 }
1025    
1026 jgs 106 double
1027 jgs 117 Data::Linf() const
1028     {
1029 jgs 123 profData->reduction1++;
1030 jgs 117 //
1031     // set the initial absolute minimum value to max double
1032 jgs 147 AbsMin abs_min_func;
1033     return algorithm(abs_min_func,numeric_limits<double>::max());
1034 jgs 117 }
1035    
1036     double
1037 jgs 106 Data::sup() const
1038 jgs 102 {
1039 jgs 123 profData->reduction1++;
1040 jgs 106 //
1041     // set the initial maximum value to min possible double
1042 jgs 147 FMax fmax_func;
1043     return algorithm(fmax_func,numeric_limits<double>::max()*-1);
1044 jgs 102 }
1045    
1046 jgs 106 double
1047     Data::inf() const
1048 jgs 102 {
1049 jgs 123 profData->reduction1++;
1050 jgs 106 //
1051     // set the initial minimum value to max possible double
1052 jgs 147 FMin fmin_func;
1053     return algorithm(fmin_func,numeric_limits<double>::max());
1054 jgs 102 }
1055    
1056     Data
1057 jgs 106 Data::maxval() const
1058 jgs 102 {
1059 jgs 123 profData->reduction2++;
1060 jgs 113 //
1061     // set the initial maximum value to min possible double
1062 jgs 147 FMax fmax_func;
1063     return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1064 jgs 102 }
1065    
1066     Data
1067 jgs 106 Data::minval() const
1068 jgs 102 {
1069 jgs 123 profData->reduction2++;
1070 jgs 113 //
1071     // set the initial minimum value to max possible double
1072 jgs 147 FMin fmin_func;
1073     return dp_algorithm(fmin_func,numeric_limits<double>::max());
1074 jgs 102 }
1075    
1076 jgs 123 Data
1077     Data::length() const
1078     {
1079     profData->reduction2++;
1080 jgs 147 Length len_func;
1081     return dp_algorithm(len_func,0);
1082 jgs 123 }
1083    
1084     Data
1085     Data::trace() const
1086     {
1087     profData->reduction2++;
1088 jgs 147 Trace trace_func;
1089     return dp_algorithm(trace_func,0);
1090 jgs 123 }
1091    
1092     Data
1093     Data::transpose(int axis) const
1094     {
1095     profData->reduction2++;
1096     // not implemented
1097     throw DataException("Error - Data::transpose not implemented yet.");
1098     return Data();
1099     }
1100    
1101 jgs 121 const boost::python::tuple
1102     Data::mindp() const
1103     {
1104 jgs 148 // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1105     // abort (for unknown reasons) if there are openmp directives with it in the
1106     // surrounding function
1107    
1108     int SampleNo;
1109     int DataPointNo;
1110    
1111     calc_mindp(SampleNo,DataPointNo);
1112    
1113     return make_tuple(SampleNo,DataPointNo);
1114     }
1115    
1116     void
1117     Data::calc_mindp(int& SampleNo,
1118     int& DataPointNo) const
1119     {
1120     int i,j;
1121     int lowi=0,lowj=0;
1122     double min=numeric_limits<double>::max();
1123    
1124 jgs 121 Data temp=minval();
1125    
1126     int numSamples=temp.getNumSamples();
1127     int numDPPSample=temp.getNumDataPointsPerSample();
1128    
1129 jgs 148 double next,local_min;
1130     int local_lowi,local_lowj;
1131 jgs 121
1132 jgs 148 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1133     {
1134     local_min=min;
1135     #pragma omp for private(i,j) schedule(static)
1136     for (i=0; i<numSamples; i++) {
1137     for (j=0; j<numDPPSample; j++) {
1138     next=temp.getDataPoint(i,j)();
1139     if (next<local_min) {
1140     local_min=next;
1141     local_lowi=i;
1142     local_lowj=j;
1143     }
1144 jgs 121 }
1145     }
1146 jgs 148 #pragma omp critical
1147     if (local_min<min) {
1148     min=local_min;
1149     lowi=local_lowi;
1150     lowj=local_lowj;
1151     }
1152 jgs 121 }
1153    
1154 jgs 148 SampleNo = lowi;
1155     DataPointNo = lowj;
1156 jgs 121 }
1157    
1158 jgs 104 void
1159     Data::saveDX(std::string fileName) const
1160     {
1161     getDomain().saveDX(fileName,*this);
1162     return;
1163     }
1164    
1165 jgs 110 void
1166     Data::saveVTK(std::string fileName) const
1167     {
1168     getDomain().saveVTK(fileName,*this);
1169     return;
1170     }
1171    
1172 jgs 102 Data&
1173     Data::operator+=(const Data& right)
1174     {
1175 jgs 123 profData->binary++;
1176 jgs 94 binaryOp(right,plus<double>());
1177     return (*this);
1178     }
1179    
1180 jgs 102 Data&
1181     Data::operator+=(const boost::python::object& right)
1182 jgs 94 {
1183 jgs 123 profData->binary++;
1184 jgs 94 binaryOp(right,plus<double>());
1185     return (*this);
1186     }
1187    
1188 jgs 102 Data&
1189     Data::operator-=(const Data& right)
1190 jgs 94 {
1191 jgs 123 profData->binary++;
1192 jgs 94 binaryOp(right,minus<double>());
1193     return (*this);
1194     }
1195    
1196 jgs 102 Data&
1197     Data::operator-=(const boost::python::object& right)
1198 jgs 94 {
1199 jgs 123 profData->binary++;
1200 jgs 94 binaryOp(right,minus<double>());
1201     return (*this);
1202     }
1203    
1204 jgs 102 Data&
1205     Data::operator*=(const Data& right)
1206 jgs 94 {
1207 jgs 123 profData->binary++;
1208 jgs 94 binaryOp(right,multiplies<double>());
1209     return (*this);
1210     }
1211    
1212 jgs 102 Data&
1213     Data::operator*=(const boost::python::object& right)
1214 jgs 94 {
1215 jgs 123 profData->binary++;
1216 jgs 94 binaryOp(right,multiplies<double>());
1217     return (*this);
1218     }
1219    
1220 jgs 102 Data&
1221     Data::operator/=(const Data& right)
1222 jgs 94 {
1223 jgs 123 profData->binary++;
1224 jgs 94 binaryOp(right,divides<double>());
1225     return (*this);
1226     }
1227    
1228 jgs 102 Data&
1229     Data::operator/=(const boost::python::object& right)
1230 jgs 94 {
1231 jgs 123 profData->binary++;
1232 jgs 94 binaryOp(right,divides<double>());
1233     return (*this);
1234     }
1235    
1236 jgs 102 Data
1237     Data::powO(const boost::python::object& right) const
1238 jgs 94 {
1239 jgs 123 profData->binary++;
1240 jgs 94 Data result;
1241     result.copy(*this);
1242     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1243     return result;
1244     }
1245    
1246 jgs 102 Data
1247     Data::powD(const Data& right) const
1248 jgs 94 {
1249 jgs 123 profData->binary++;
1250 jgs 94 Data result;
1251     result.copy(*this);
1252     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1253     return result;
1254     }
1255    
1256     //
1257 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1258 jgs 102 Data
1259     escript::operator+(const Data& left, const Data& right)
1260 jgs 94 {
1261     Data result;
1262     //
1263     // perform a deep copy
1264     result.copy(left);
1265     result+=right;
1266     return result;
1267     }
1268    
1269     //
1270 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1271 jgs 102 Data
1272     escript::operator-(const Data& left, const Data& right)
1273 jgs 94 {
1274     Data result;
1275     //
1276     // perform a deep copy
1277     result.copy(left);
1278     result-=right;
1279     return result;
1280     }
1281    
1282     //
1283 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1284 jgs 102 Data
1285     escript::operator*(const Data& left, const Data& right)
1286 jgs 94 {
1287     Data result;
1288     //
1289     // perform a deep copy
1290     result.copy(left);
1291     result*=right;
1292     return result;
1293     }
1294    
1295     //
1296 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1297 jgs 102 Data
1298     escript::operator/(const Data& left, const Data& right)
1299 jgs 94 {
1300     Data result;
1301     //
1302     // perform a deep copy
1303     result.copy(left);
1304     result/=right;
1305     return result;
1306     }
1307    
1308     //
1309 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1310 jgs 102 Data
1311     escript::operator+(const Data& left, const boost::python::object& right)
1312 jgs 94 {
1313     //
1314     // Convert to DataArray format if possible
1315     DataArray temp(right);
1316     Data result;
1317     //
1318     // perform a deep copy
1319     result.copy(left);
1320     result+=right;
1321     return result;
1322     }
1323    
1324     //
1325 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1326 jgs 102 Data
1327     escript::operator-(const Data& left, const boost::python::object& right)
1328 jgs 94 {
1329     //
1330     // Convert to DataArray format if possible
1331     DataArray temp(right);
1332     Data result;
1333     //
1334     // perform a deep copy
1335     result.copy(left);
1336     result-=right;
1337     return result;
1338     }
1339    
1340     //
1341 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1342 jgs 102 Data
1343     escript::operator*(const Data& left, const boost::python::object& right)
1344 jgs 94 {
1345     //
1346     // Convert to DataArray format if possible
1347     DataArray temp(right);
1348     Data result;
1349     //
1350     // perform a deep copy
1351     result.copy(left);
1352     result*=right;
1353     return result;
1354     }
1355    
1356     //
1357 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1358 jgs 102 Data
1359     escript::operator/(const Data& left, const boost::python::object& right)
1360 jgs 94 {
1361     //
1362     // Convert to DataArray format if possible
1363     DataArray temp(right);
1364     Data result;
1365     //
1366     // perform a deep copy
1367     result.copy(left);
1368     result/=right;
1369     return result;
1370     }
1371    
1372     //
1373 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1374 jgs 102 Data
1375     escript::operator+(const boost::python::object& left, const Data& right)
1376 jgs 94 {
1377     //
1378     // Construct the result using the given value and the other parameters
1379     // from right
1380     Data result(left,right);
1381     result+=right;
1382     return result;
1383     }
1384    
1385     //
1386 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1387 jgs 102 Data
1388     escript::operator-(const boost::python::object& left, const Data& right)
1389 jgs 94 {
1390     //
1391     // Construct the result using the given value and the other parameters
1392     // from right
1393     Data result(left,right);
1394     result-=right;
1395     return result;
1396     }
1397    
1398     //
1399 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1400 jgs 102 Data
1401     escript::operator*(const boost::python::object& left, const Data& right)
1402 jgs 94 {
1403     //
1404     // Construct the result using the given value and the other parameters
1405     // from right
1406     Data result(left,right);
1407     result*=right;
1408     return result;
1409     }
1410    
1411     //
1412 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1413 jgs 102 Data
1414     escript::operator/(const boost::python::object& left, const Data& right)
1415 jgs 94 {
1416     //
1417     // Construct the result using the given value and the other parameters
1418     // from right
1419     Data result(left,right);
1420     result/=right;
1421     return result;
1422     }
1423    
1424     //
1425 jgs 102 //bool escript::operator==(const Data& left, const Data& right)
1426     //{
1427     // /*
1428     // NB: this operator does very little at this point, and isn't to
1429     // be relied on. Requires further implementation.
1430     // */
1431     //
1432     // bool ret;
1433     //
1434     // if (left.isEmpty()) {
1435     // if(!right.isEmpty()) {
1436     // ret = false;
1437     // } else {
1438     // ret = true;
1439     // }
1440     // }
1441     //
1442     // if (left.isConstant()) {
1443     // if(!right.isConstant()) {
1444     // ret = false;
1445     // } else {
1446     // ret = true;
1447     // }
1448     // }
1449     //
1450     // if (left.isTagged()) {
1451     // if(!right.isTagged()) {
1452     // ret = false;
1453     // } else {
1454     // ret = true;
1455     // }
1456     // }
1457     //
1458     // if (left.isExpanded()) {
1459     // if(!right.isExpanded()) {
1460     // ret = false;
1461     // } else {
1462     // ret = true;
1463     // }
1464     // }
1465     //
1466     // return ret;
1467     //}
1468    
1469     Data
1470     Data::getItem(const boost::python::object& key) const
1471 jgs 94 {
1472 jgs 102 const DataArrayView& view=getPointDataView();
1473 jgs 94
1474 jgs 102 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1475 jgs 94
1476 jgs 102 if (slice_region.size()!=view.getRank()) {
1477     throw DataException("Error - slice size does not match Data rank.");
1478 jgs 94 }
1479    
1480 jgs 102 return getSlice(slice_region);
1481 jgs 94 }
1482    
1483     Data
1484 jgs 102 Data::getSlice(const DataArrayView::RegionType& region) const
1485 jgs 94 {
1486 jgs 123 profData->slicing++;
1487 jgs 102 return Data(*this,region);
1488 jgs 94 }
1489    
1490     void
1491 jgs 102 Data::setItemO(const boost::python::object& key,
1492     const boost::python::object& value)
1493 jgs 94 {
1494 jgs 102 Data tempData(value,getFunctionSpace());
1495     setItemD(key,tempData);
1496     }
1497    
1498     void
1499     Data::setItemD(const boost::python::object& key,
1500     const Data& value)
1501     {
1502 jgs 94 const DataArrayView& view=getPointDataView();
1503 jgs 123
1504 jgs 94 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1505     if (slice_region.size()!=view.getRank()) {
1506     throw DataException("Error - slice size does not match Data rank.");
1507     }
1508 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
1509     setSlice(Data(value,getFunctionSpace()),slice_region);
1510     } else {
1511     setSlice(value,slice_region);
1512     }
1513 jgs 94 }
1514    
1515     void
1516 jgs 102 Data::setSlice(const Data& value,
1517     const DataArrayView::RegionType& region)
1518 jgs 94 {
1519 jgs 123 profData->slicing++;
1520 jgs 102 Data tempValue(value);
1521     typeMatchLeft(tempValue);
1522     typeMatchRight(tempValue);
1523     m_data->setSlice(tempValue.m_data.get(),region);
1524     }
1525    
1526     void
1527     Data::typeMatchLeft(Data& right) const
1528     {
1529     if (isExpanded()){
1530     right.expand();
1531     } else if (isTagged()) {
1532     if (right.isConstant()) {
1533     right.tag();
1534     }
1535     }
1536     }
1537    
1538     void
1539     Data::typeMatchRight(const Data& right)
1540     {
1541 jgs 94 if (isTagged()) {
1542     if (right.isExpanded()) {
1543     expand();
1544     }
1545     } else if (isConstant()) {
1546     if (right.isExpanded()) {
1547     expand();
1548     } else if (right.isTagged()) {
1549     tag();
1550     }
1551     }
1552     }
1553    
1554     void
1555     Data::setTaggedValue(int tagKey,
1556     const boost::python::object& value)
1557     {
1558     //
1559     // Ensure underlying data object is of type DataTagged
1560     tag();
1561    
1562     if (!isTagged()) {
1563     throw DataException("Error - DataTagged conversion failed!!");
1564     }
1565    
1566     //
1567     // Construct DataArray from boost::python::object input value
1568     DataArray valueDataArray(value);
1569    
1570     //
1571     // Call DataAbstract::setTaggedValue
1572     m_data->setTaggedValue(tagKey,valueDataArray.getView());
1573     }
1574    
1575 jgs 110 void
1576 jgs 121 Data::setTaggedValueFromCPP(int tagKey,
1577     const DataArrayView& value)
1578     {
1579     //
1580     // Ensure underlying data object is of type DataTagged
1581     tag();
1582    
1583     if (!isTagged()) {
1584     throw DataException("Error - DataTagged conversion failed!!");
1585     }
1586    
1587     //
1588     // Call DataAbstract::setTaggedValue
1589     m_data->setTaggedValue(tagKey,value);
1590     }
1591    
1592 jgs 149 int
1593     Data::getTagNumber(int dpno)
1594     {
1595     return m_data->getTagNumber(dpno);
1596     }
1597    
1598 jgs 121 void
1599 jgs 110 Data::setRefValue(int ref,
1600     const boost::python::numeric::array& value)
1601     {
1602     //
1603     // Construct DataArray from boost::python::object input value
1604     DataArray valueDataArray(value);
1605    
1606     //
1607     // Call DataAbstract::setRefValue
1608     m_data->setRefValue(ref,valueDataArray);
1609     }
1610    
1611     void
1612     Data::getRefValue(int ref,
1613     boost::python::numeric::array& value)
1614     {
1615     //
1616     // Construct DataArray for boost::python::object return value
1617     DataArray valueDataArray(value);
1618    
1619     //
1620     // Load DataArray with values from data-points specified by ref
1621     m_data->getRefValue(ref,valueDataArray);
1622    
1623     //
1624     // Load values from valueDataArray into return numarray
1625    
1626     // extract the shape of the numarray
1627     int rank = value.getrank();
1628     DataArrayView::ShapeType shape;
1629     for (int i=0; i < rank; i++) {
1630     shape.push_back(extract<int>(value.getshape()[i]));
1631     }
1632    
1633     // and load the numarray with the data from the DataArray
1634     DataArrayView valueView = valueDataArray.getView();
1635    
1636     if (rank==0) {
1637 jgs 126 boost::python::numeric::array temp_numArray(valueView());
1638     value = temp_numArray;
1639 jgs 110 }
1640     if (rank==1) {
1641     for (int i=0; i < shape[0]; i++) {
1642     value[i] = valueView(i);
1643     }
1644     }
1645     if (rank==2) {
1646 jgs 126 for (int i=0; i < shape[0]; i++) {
1647     for (int j=0; j < shape[1]; j++) {
1648     value[i][j] = valueView(i,j);
1649     }
1650     }
1651 jgs 110 }
1652     if (rank==3) {
1653 jgs 126 for (int i=0; i < shape[0]; i++) {
1654     for (int j=0; j < shape[1]; j++) {
1655     for (int k=0; k < shape[2]; k++) {
1656     value[i][j][k] = valueView(i,j,k);
1657     }
1658     }
1659     }
1660 jgs 110 }
1661     if (rank==4) {
1662 jgs 126 for (int i=0; i < shape[0]; i++) {
1663     for (int j=0; j < shape[1]; j++) {
1664     for (int k=0; k < shape[2]; k++) {
1665     for (int l=0; l < shape[3]; l++) {
1666     value[i][j][k][l] = valueView(i,j,k,l);
1667     }
1668     }
1669     }
1670     }
1671 jgs 110 }
1672    
1673     }
1674    
1675 jgs 94 void
1676 jgs 119 Data::archiveData(const std::string fileName)
1677     {
1678     cout << "Archiving Data object to: " << fileName << endl;
1679    
1680     //
1681     // Determine type of this Data object
1682     int dataType = -1;
1683    
1684     if (isEmpty()) {
1685     dataType = 0;
1686     cout << "\tdataType: DataEmpty" << endl;
1687     }
1688     if (isConstant()) {
1689     dataType = 1;
1690     cout << "\tdataType: DataConstant" << endl;
1691     }
1692     if (isTagged()) {
1693     dataType = 2;
1694     cout << "\tdataType: DataTagged" << endl;
1695     }
1696     if (isExpanded()) {
1697     dataType = 3;
1698     cout << "\tdataType: DataExpanded" << endl;
1699     }
1700 jgs 123
1701 jgs 119 if (dataType == -1) {
1702     throw DataException("archiveData Error: undefined dataType");
1703     }
1704    
1705     //
1706     // Collect data items common to all Data types
1707     int noSamples = getNumSamples();
1708     int noDPPSample = getNumDataPointsPerSample();
1709     int functionSpaceType = getFunctionSpace().getTypeCode();
1710     int dataPointRank = getDataPointRank();
1711     int dataPointSize = getDataPointSize();
1712     int dataLength = getLength();
1713     DataArrayView::ShapeType dataPointShape = getDataPointShape();
1714     int referenceNumbers[noSamples];
1715     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1716     referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1717     }
1718     int tagNumbers[noSamples];
1719     if (isTagged()) {
1720     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1721     tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1722     }
1723     }
1724    
1725     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1726     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1727     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1728    
1729     //
1730     // Flatten Shape to an array of integers suitable for writing to file
1731     int flatShape[4] = {0,0,0,0};
1732     cout << "\tshape: < ";
1733     for (int dim=0; dim<dataPointRank; dim++) {
1734     flatShape[dim] = dataPointShape[dim];
1735     cout << dataPointShape[dim] << " ";
1736     }
1737     cout << ">" << endl;
1738    
1739     //
1740 jgs 123 // Open archive file
1741 jgs 119 ofstream archiveFile;
1742     archiveFile.open(fileName.data(), ios::out);
1743    
1744     if (!archiveFile.good()) {
1745     throw DataException("archiveData Error: problem opening archive file");
1746     }
1747    
1748 jgs 123 //
1749     // Write common data items to archive file
1750 jgs 119 archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1751     archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1752     archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1753     archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1754     archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1755     archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1756     archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1757     for (int dim = 0; dim < 4; dim++) {
1758     archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1759     }
1760     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1761     archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1762     }
1763     if (isTagged()) {
1764     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1765     archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1766     }
1767     }
1768    
1769     if (!archiveFile.good()) {
1770     throw DataException("archiveData Error: problem writing to archive file");
1771     }
1772    
1773     //
1774 jgs 123 // Archive underlying data values for each Data type
1775     int noValues;
1776 jgs 119 switch (dataType) {
1777     case 0:
1778     // DataEmpty
1779 jgs 123 noValues = 0;
1780     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1781     cout << "\tnoValues: " << noValues << endl;
1782 jgs 119 break;
1783     case 1:
1784     // DataConstant
1785 jgs 123 noValues = m_data->getLength();
1786     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1787     cout << "\tnoValues: " << noValues << endl;
1788     if (m_data->archiveData(archiveFile,noValues)) {
1789     throw DataException("archiveData Error: problem writing data to archive file");
1790     }
1791 jgs 119 break;
1792     case 2:
1793     // DataTagged
1794 jgs 123 noValues = m_data->getLength();
1795     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1796     cout << "\tnoValues: " << noValues << endl;
1797     if (m_data->archiveData(archiveFile,noValues)) {
1798     throw DataException("archiveData Error: problem writing data to archive file");
1799     }
1800 jgs 119 break;
1801     case 3:
1802     // DataExpanded
1803 jgs 123 noValues = m_data->getLength();
1804     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1805     cout << "\tnoValues: " << noValues << endl;
1806     if (m_data->archiveData(archiveFile,noValues)) {
1807     throw DataException("archiveData Error: problem writing data to archive file");
1808     }
1809 jgs 119 break;
1810     }
1811    
1812 jgs 123 if (!archiveFile.good()) {
1813     throw DataException("archiveData Error: problem writing data to archive file");
1814     }
1815    
1816     //
1817     // Close archive file
1818     archiveFile.close();
1819    
1820     if (!archiveFile.good()) {
1821     throw DataException("archiveData Error: problem closing archive file");
1822     }
1823    
1824 jgs 119 }
1825    
1826     void
1827     Data::extractData(const std::string fileName,
1828     const FunctionSpace& fspace)
1829     {
1830     //
1831     // Can only extract Data to an object which is initially DataEmpty
1832     if (!isEmpty()) {
1833     throw DataException("extractData Error: can only extract to DataEmpty object");
1834     }
1835    
1836     cout << "Extracting Data object from: " << fileName << endl;
1837    
1838     int dataType;
1839     int noSamples;
1840     int noDPPSample;
1841     int functionSpaceType;
1842     int dataPointRank;
1843     int dataPointSize;
1844     int dataLength;
1845     DataArrayView::ShapeType dataPointShape;
1846     int flatShape[4];
1847    
1848     //
1849 jgs 123 // Open the archive file
1850 jgs 119 ifstream archiveFile;
1851     archiveFile.open(fileName.data(), ios::in);
1852    
1853     if (!archiveFile.good()) {
1854     throw DataException("extractData Error: problem opening archive file");
1855     }
1856    
1857 jgs 123 //
1858     // Read common data items from archive file
1859 jgs 119 archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1860     archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1861     archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1862     archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1863     archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1864     archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1865     archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1866     for (int dim = 0; dim < 4; dim++) {
1867     archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1868     if (flatShape[dim]>0) {
1869     dataPointShape.push_back(flatShape[dim]);
1870     }
1871     }
1872     int referenceNumbers[noSamples];
1873     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1874     archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1875     }
1876     int tagNumbers[noSamples];
1877     if (dataType==2) {
1878     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1879     archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1880     }
1881     }
1882    
1883     if (!archiveFile.good()) {
1884     throw DataException("extractData Error: problem reading from archive file");
1885     }
1886    
1887 jgs 123 //
1888     // Verify the values just read from the archive file
1889 jgs 119 switch (dataType) {
1890     case 0:
1891     cout << "\tdataType: DataEmpty" << endl;
1892     break;
1893     case 1:
1894     cout << "\tdataType: DataConstant" << endl;
1895     break;
1896     case 2:
1897     cout << "\tdataType: DataTagged" << endl;
1898     break;
1899     case 3:
1900     cout << "\tdataType: DataExpanded" << endl;
1901     break;
1902     default:
1903     throw DataException("extractData Error: undefined dataType read from archive file");
1904     break;
1905     }
1906    
1907     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1908     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1909     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1910     cout << "\tshape: < ";
1911     for (int dim = 0; dim < dataPointRank; dim++) {
1912     cout << dataPointShape[dim] << " ";
1913     }
1914     cout << ">" << endl;
1915    
1916     //
1917     // Verify that supplied FunctionSpace object is compatible with this Data object.
1918     if ( (fspace.getTypeCode()!=functionSpaceType) ||
1919     (fspace.getNumSamples()!=noSamples) ||
1920     (fspace.getNumDPPSample()!=noDPPSample)
1921     ) {
1922     throw DataException("extractData Error: incompatible FunctionSpace");
1923     }
1924     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1925     if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
1926     throw DataException("extractData Error: incompatible FunctionSpace");
1927     }
1928     }
1929     if (dataType==2) {
1930     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1931     if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
1932     throw DataException("extractData Error: incompatible FunctionSpace");
1933     }
1934     }
1935     }
1936    
1937     //
1938     // Construct a DataVector to hold underlying data values
1939     DataVector dataVec(dataLength);
1940    
1941     //
1942     // Load this DataVector with the appropriate values
1943 jgs 123 int noValues;
1944     archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
1945     cout << "\tnoValues: " << noValues << endl;
1946 jgs 119 switch (dataType) {
1947     case 0:
1948     // DataEmpty
1949 jgs 123 if (noValues != 0) {
1950     throw DataException("extractData Error: problem reading data from archive file");
1951     }
1952 jgs 119 break;
1953     case 1:
1954     // DataConstant
1955 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
1956     throw DataException("extractData Error: problem reading data from archive file");
1957     }
1958 jgs 119 break;
1959     case 2:
1960     // DataTagged
1961 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
1962     throw DataException("extractData Error: problem reading data from archive file");
1963     }
1964 jgs 119 break;
1965     case 3:
1966     // DataExpanded
1967 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
1968     throw DataException("extractData Error: problem reading data from archive file");
1969     }
1970 jgs 119 break;
1971     }
1972    
1973 jgs 123 if (!archiveFile.good()) {
1974     throw DataException("extractData Error: problem reading from archive file");
1975     }
1976    
1977 jgs 119 //
1978 jgs 123 // Close archive file
1979     archiveFile.close();
1980    
1981     if (!archiveFile.good()) {
1982     throw DataException("extractData Error: problem closing archive file");
1983     }
1984    
1985     //
1986 jgs 119 // Construct an appropriate Data object
1987     DataAbstract* tempData;
1988     switch (dataType) {
1989     case 0:
1990     // DataEmpty
1991     tempData=new DataEmpty();
1992     break;
1993     case 1:
1994     // DataConstant
1995     tempData=new DataConstant(fspace,dataPointShape,dataVec);
1996     break;
1997     case 2:
1998     // DataTagged
1999     tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2000     break;
2001     case 3:
2002     // DataExpanded
2003     tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2004     break;
2005     }
2006     shared_ptr<DataAbstract> temp_data(tempData);
2007     m_data=temp_data;
2008     }
2009    
2010 jgs 94 ostream& escript::operator<<(ostream& o, const Data& data)
2011     {
2012     o << data.toString();
2013     return o;
2014     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26