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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26