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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 119 - (hide annotations)
Tue Apr 12 04:45:05 2005 UTC (14 years, 7 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/Data.cpp
File size: 37286 byte(s)
*** empty log message ***

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     #include "escript/Data/DataAlgorithm.h"
39     #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     Data::Data()
49     {
50     //
51     // Default data is type DataEmpty
52     DataAbstract* temp=new DataEmpty();
53 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
54     m_data=temp_data;
55 jgs 94 }
56    
57     Data::Data(double value,
58     const tuple& shape,
59     const FunctionSpace& what,
60     bool expanded)
61     {
62     DataArrayView::ShapeType dataPointShape;
63     for (int i = 0; i < shape.attr("__len__")(); ++i) {
64     dataPointShape.push_back(extract<const int>(shape[i]));
65     }
66     DataArray temp(dataPointShape,value);
67     initialise(temp.getView(),what,expanded);
68     }
69    
70     Data::Data(double value,
71     const DataArrayView::ShapeType& dataPointShape,
72     const FunctionSpace& what,
73     bool expanded)
74     {
75     DataArray temp(dataPointShape,value);
76     pair<int,int> dataShape=what.getDataShape();
77     initialise(temp.getView(),what,expanded);
78     }
79    
80 jgs 102 Data::Data(const Data& inData)
81 jgs 94 {
82 jgs 102 m_data=inData.m_data;
83 jgs 94 }
84    
85     Data::Data(const Data& inData,
86     const DataArrayView::RegionType& region)
87     {
88     //
89 jgs 102 // Create Data which is a slice of another Data
90     DataAbstract* tmp = inData.m_data->getSlice(region);
91     shared_ptr<DataAbstract> temp_data(tmp);
92     m_data=temp_data;
93 jgs 94 }
94    
95     Data::Data(const Data& inData,
96     const FunctionSpace& functionspace)
97     {
98     if (inData.getFunctionSpace()==functionspace) {
99     m_data=inData.m_data;
100     } else {
101     Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
102     // Note for Lutz, Must use a reference or pointer to a derived object
103     // in order to get polymorphic behaviour. Shouldn't really
104     // be able to create an instance of AbstractDomain but that was done
105     // as a boost python work around which may no longer be required.
106     const AbstractDomain& inDataDomain=inData.getDomain();
107     if (inDataDomain==functionspace.getDomain()) {
108     inDataDomain.interpolateOnDomain(tmp,inData);
109     } else {
110     inDataDomain.interpolateACross(tmp,inData);
111     }
112     m_data=tmp.m_data;
113     }
114     }
115    
116     Data::Data(const DataTagged::TagListType& tagKeys,
117     const DataTagged::ValueListType & values,
118     const DataArrayView& defaultValue,
119     const FunctionSpace& what,
120     bool expanded)
121     {
122     DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
123 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
124     m_data=temp_data;
125 jgs 94 if (expanded) {
126     expand();
127     }
128     }
129    
130     Data::Data(const numeric::array& value,
131     const FunctionSpace& what,
132     bool expanded)
133     {
134     initialise(value,what,expanded);
135     }
136    
137     Data::Data(const DataArrayView& value,
138     const FunctionSpace& what,
139     bool expanded)
140     {
141     initialise(value,what,expanded);
142     }
143    
144     Data::Data(const object& value,
145     const FunctionSpace& what,
146     bool expanded)
147     {
148     numeric::array asNumArray(value);
149     initialise(asNumArray,what,expanded);
150     }
151    
152     Data::Data(const object& value,
153     const Data& other)
154     {
155     //
156     // Create DataConstant using the given value and all other parameters
157     // copied from other. If value is a rank 0 object this Data
158     // will assume the point data shape of other.
159     DataArray temp(value);
160     if (temp.getView().getRank()==0) {
161     //
162     // Create a DataArray with the scalar value for all elements
163     DataArray temp2(other.getPointDataView().getShape(),temp.getView()());
164     initialise(temp2.getView(),other.getFunctionSpace(),false);
165     } else {
166     //
167     // Create a DataConstant with the same sample shape as other
168     initialise(temp.getView(),other.getFunctionSpace(),false);
169     }
170     }
171    
172     escriptDataC
173     Data::getDataC()
174     {
175     escriptDataC temp;
176     temp.m_dataPtr=(void*)this;
177     return temp;
178     }
179    
180     escriptDataC
181     Data::getDataC() const
182     {
183     escriptDataC temp;
184     temp.m_dataPtr=(void*)this;
185     return temp;
186     }
187    
188     tuple
189     Data::getShapeTuple() const
190     {
191     const DataArrayView::ShapeType& shape=getDataPointShape();
192     switch(getDataPointRank()) {
193     case 0:
194     return make_tuple();
195     case 1:
196     return make_tuple(long_(shape[0]));
197     case 2:
198     return make_tuple(long_(shape[0]),long_(shape[1]));
199     case 3:
200     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]));
201     case 4:
202     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]),long_(shape[3]));
203     default:
204     throw DataException("Error - illegal Data rank.");
205     }
206     }
207    
208     void
209     Data::copy(const Data& other)
210     {
211     //
212     // Perform a deep copy
213     {
214     DataExpanded* temp=dynamic_cast<DataExpanded*>(other.m_data.get());
215     if (temp!=0) {
216     //
217     // Construct a DataExpanded copy
218     DataAbstract* newData=new DataExpanded(*temp);
219 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
220     m_data=temp_data;
221 jgs 94 return;
222     }
223     }
224     {
225     DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());
226     if (temp!=0) {
227     //
228 jgs 102 // Construct a DataTagged copy
229 jgs 94 DataAbstract* newData=new DataTagged(*temp);
230 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
231     m_data=temp_data;
232 jgs 94 return;
233     }
234     }
235     {
236     DataConstant* temp=dynamic_cast<DataConstant*>(other.m_data.get());
237     if (temp!=0) {
238     //
239     // Construct a DataConstant copy
240     DataAbstract* newData=new DataConstant(*temp);
241 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
242     m_data=temp_data;
243 jgs 94 return;
244     }
245     }
246 jgs 102 {
247     DataEmpty* temp=dynamic_cast<DataEmpty*>(other.m_data.get());
248     if (temp!=0) {
249     //
250     // Construct a DataEmpty copy
251     DataAbstract* newData=new DataEmpty();
252     shared_ptr<DataAbstract> temp_data(newData);
253     m_data=temp_data;
254     return;
255     }
256     }
257 jgs 94 throw DataException("Error - Copy not implemented for this Data type.");
258     }
259    
260     void
261     Data::copyWithMask(const Data& other,
262     const Data& mask)
263     {
264     Data mask1;
265     Data mask2;
266    
267     mask1 = mask.wherePositive();
268     mask2.copy(mask1);
269    
270     mask1 *= other;
271     mask2 *= *this;
272     mask2 = *this - mask2;
273    
274     *this = mask1 + mask2;
275     }
276    
277     bool
278     Data::isExpanded() const
279     {
280     DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
281     return (temp!=0);
282     }
283    
284     bool
285     Data::isTagged() const
286     {
287     DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
288     return (temp!=0);
289     }
290    
291     bool
292     Data::isEmpty() const
293     {
294     DataEmpty* temp=dynamic_cast<DataEmpty*>(m_data.get());
295     return (temp!=0);
296     }
297    
298     bool
299     Data::isConstant() const
300     {
301     DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
302     return (temp!=0);
303     }
304    
305     void
306     Data::expand()
307     {
308     if (isConstant()) {
309     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
310     DataAbstract* temp=new DataExpanded(*tempDataConst);
311 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
312     m_data=temp_data;
313 jgs 94 } else if (isTagged()) {
314     DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
315     DataAbstract* temp=new DataExpanded(*tempDataTag);
316 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
317     m_data=temp_data;
318 jgs 94 } else if (isExpanded()) {
319     //
320     // do nothing
321     } else if (isEmpty()) {
322     throw DataException("Error - Expansion of DataEmpty not possible.");
323     } else {
324     throw DataException("Error - Expansion not implemented for this Data type.");
325     }
326     }
327    
328     void
329     Data::tag()
330     {
331     if (isConstant()) {
332     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
333     DataAbstract* temp=new DataTagged(*tempDataConst);
334 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
335     m_data=temp_data;
336 jgs 94 } else if (isTagged()) {
337     // do nothing
338     } else if (isExpanded()) {
339     throw DataException("Error - Creating tag data from DataExpanded not possible.");
340     } else if (isEmpty()) {
341     throw DataException("Error - Creating tag data from DataEmpty not possible.");
342     } else {
343     throw DataException("Error - Tagging not implemented for this Data type.");
344     }
345     }
346    
347 jgs 102 void
348     Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)
349     {
350     m_data->reshapeDataPoint(shape);
351     }
352    
353 jgs 94 Data
354     Data::wherePositive() const
355     {
356     return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
357     }
358    
359     Data
360 jgs 102 Data::whereNegative() const
361     {
362     return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
363     }
364    
365     Data
366 jgs 94 Data::whereNonNegative() const
367     {
368     return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
369     }
370    
371     Data
372 jgs 102 Data::whereNonPositive() const
373 jgs 94 {
374 jgs 102 return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
375 jgs 94 }
376    
377     Data
378     Data::whereZero() const
379     {
380     return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));
381     }
382    
383     Data
384 jgs 102 Data::whereNonZero() const
385     {
386     return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));
387     }
388    
389     Data
390 jgs 94 Data::interpolate(const FunctionSpace& functionspace) const
391     {
392     return Data(*this,functionspace);
393     }
394    
395     bool
396     Data::probeInterpolation(const FunctionSpace& functionspace) const
397     {
398     if (getFunctionSpace()==functionspace) {
399     return true;
400     } else {
401     const AbstractDomain& domain=getDomain();
402     if (domain==functionspace.getDomain()) {
403     return domain.probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());
404     } else {
405     return domain.probeInterpolationACross(getFunctionSpace().getTypeCode(),functionspace.getDomain(),functionspace.getTypeCode());
406     }
407     }
408     }
409    
410     Data
411     Data::gradOn(const FunctionSpace& functionspace) const
412     {
413     if (functionspace.getDomain()!=getDomain())
414     throw DataException("Error - gradient cannot be calculated on different domains.");
415     DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
416     grad_shape.push_back(functionspace.getDim());
417     Data out(0.0,grad_shape,functionspace,true);
418     getDomain().setToGradient(out,*this);
419     return out;
420     }
421    
422     Data
423     Data::grad() const
424     {
425     return gradOn(escript::function(getDomain()));
426     }
427    
428     int
429     Data::getDataPointSize() const
430     {
431     return getPointDataView().noValues();
432     }
433    
434     DataArrayView::ValueType::size_type
435     Data::getLength() const
436     {
437     return m_data->getLength();
438     }
439    
440     const DataArrayView::ShapeType&
441     Data::getDataPointShape() const
442     {
443     return getPointDataView().getShape();
444     }
445    
446     boost::python::numeric::array
447 jgs 117 Data::convertToNumArray()
448     {
449     //
450     // determine the current number of data points
451     int numSamples = getNumSamples();
452     int numDataPointsPerSample = getNumDataPointsPerSample();
453     int numDataPoints = numSamples * numDataPointsPerSample;
454    
455     //
456     // determine the rank and shape of each data point
457     int dataPointRank = getDataPointRank();
458     DataArrayView::ShapeType dataPointShape = getDataPointShape();
459    
460     //
461     // create the numeric array to be returned
462     boost::python::numeric::array numArray(0.0);
463    
464     //
465     // the rank of the returned numeric array will be the rank of
466     // the data points, plus one. Where the rank of the array is n,
467     // the last n-1 dimensions will be equal to the shape of the
468     // data points, whilst the first dimension will be equal to the
469     // total number of data points. Thus the array will consist of
470     // a serial vector of the data points.
471     int arrayRank = dataPointRank + 1;
472     DataArrayView::ShapeType arrayShape;
473     arrayShape.push_back(numDataPoints);
474     for (int d=0; d<dataPointRank; d++) {
475     arrayShape.push_back(dataPointShape[d]);
476     }
477    
478     //
479     // resize the numeric array to the shape just calculated
480     if (arrayRank==1) {
481     numArray.resize(arrayShape[0]);
482     }
483     if (arrayRank==2) {
484     numArray.resize(arrayShape[0],arrayShape[1]);
485     }
486     if (arrayRank==3) {
487     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
488     }
489     if (arrayRank==4) {
490     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
491     }
492     if (arrayRank==5) {
493     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
494     }
495    
496     //
497     // loop through each data point in turn, loading the values for that data point
498     // into the numeric array.
499     int dataPoint = 0;
500     for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
501     for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
502    
503     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
504    
505     if (dataPointRank==0) {
506     numArray[dataPoint]=dataPointView();
507     }
508    
509     if (dataPointRank==1) {
510     for (int i=0; i<dataPointShape[0]; i++) {
511     numArray[dataPoint][i]=dataPointView(i);
512     }
513     }
514    
515     if (dataPointRank==2) {
516     for (int i=0; i<dataPointShape[0]; i++) {
517     for (int j=0; j<dataPointShape[1]; j++) {
518     numArray[dataPoint][i][j] = dataPointView(i,j);
519     }
520     }
521     }
522    
523     if (dataPointRank==3) {
524     for (int i=0; i<dataPointShape[0]; i++) {
525     for (int j=0; j<dataPointShape[1]; j++) {
526     for (int k=0; k<dataPointShape[2]; k++) {
527     numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
528     }
529     }
530     }
531     }
532    
533     if (dataPointRank==4) {
534     for (int i=0; i<dataPointShape[0]; i++) {
535     for (int j=0; j<dataPointShape[1]; j++) {
536     for (int k=0; k<dataPointShape[2]; k++) {
537     for (int l=0; l<dataPointShape[3]; l++) {
538     numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
539     }
540     }
541     }
542     }
543     }
544    
545     dataPoint++;
546    
547     }
548     }
549    
550     //
551     // return the loaded array
552     return numArray;
553     }
554    
555     boost::python::numeric::array
556 jgs 94 Data::integrate() const
557     {
558     int index;
559     int rank = getDataPointRank();
560     DataArrayView::ShapeType shape = getDataPointShape();
561    
562     //
563     // calculate the integral values
564     vector<double> integrals(getDataPointSize());
565     AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);
566    
567     //
568     // create the numeric array to be returned
569     // and load the array with the integral values
570     boost::python::numeric::array bp_array(1.0);
571     if (rank==0) {
572 jgs 108 bp_array.resize(1);
573 jgs 94 index = 0;
574     bp_array[0] = integrals[index];
575     }
576     if (rank==1) {
577     bp_array.resize(shape[0]);
578     for (int i=0; i<shape[0]; i++) {
579     index = i;
580     bp_array[i] = integrals[index];
581     }
582     }
583     if (rank==2) {
584     bp_array.resize(shape[0],shape[1]);
585     for (int i=0; i<shape[0]; i++) {
586     for (int j=0; j<shape[1]; j++) {
587     index = i + shape[0] * j;
588     bp_array[i,j] = integrals[index];
589     }
590     }
591     }
592     if (rank==3) {
593     bp_array.resize(shape[0],shape[1],shape[2]);
594     for (int i=0; i<shape[0]; i++) {
595     for (int j=0; j<shape[1]; j++) {
596     for (int k=0; k<shape[2]; k++) {
597     index = i + shape[0] * ( j + shape[1] * k );
598     bp_array[i,j,k] = integrals[index];
599     }
600     }
601     }
602     }
603     if (rank==4) {
604     bp_array.resize(shape[0],shape[1],shape[2],shape[3]);
605     for (int i=0; i<shape[0]; i++) {
606     for (int j=0; j<shape[1]; j++) {
607     for (int k=0; k<shape[2]; k++) {
608     for (int l=0; l<shape[3]; l++) {
609     index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
610     bp_array[i,j,k,l] = integrals[index];
611     }
612     }
613     }
614     }
615     }
616    
617     //
618     // return the loaded array
619     return bp_array;
620     }
621    
622     Data
623     Data::sin() const
624     {
625     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
626     }
627    
628     Data
629     Data::cos() const
630     {
631     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
632     }
633    
634     Data
635     Data::tan() const
636     {
637     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
638     }
639    
640     Data
641     Data::log() const
642     {
643     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
644     }
645    
646     Data
647     Data::ln() const
648     {
649     return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
650     }
651    
652 jgs 106 Data
653     Data::sign() const
654 jgs 94 {
655 jgs 106 return escript::unaryOp(*this,escript::fsign);
656 jgs 94 }
657    
658 jgs 106 Data
659     Data::abs() const
660 jgs 94 {
661 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
662 jgs 94 }
663    
664 jgs 106 Data
665     Data::neg() const
666 jgs 94 {
667 jgs 106 return escript::unaryOp(*this,negate<double>());
668 jgs 94 }
669    
670 jgs 102 Data
671 jgs 106 Data::pos() const
672 jgs 94 {
673 jgs 102 return (*this);
674     }
675    
676     Data
677 jgs 106 Data::exp() const
678 jgs 102 {
679 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
680 jgs 102 }
681    
682     Data
683 jgs 106 Data::sqrt() const
684 jgs 102 {
685 jgs 106 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
686 jgs 102 }
687    
688 jgs 106 double
689     Data::Lsup() const
690 jgs 102 {
691 jgs 106 //
692     // set the initial absolute maximum value to zero
693     return algorithm(DataAlgorithmAdapter<AbsMax>(0));
694 jgs 102 }
695    
696 jgs 106 double
697 jgs 117 Data::Linf() const
698     {
699     //
700     // set the initial absolute minimum value to max double
701     return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));
702     }
703    
704     double
705 jgs 106 Data::sup() const
706 jgs 102 {
707 jgs 106 //
708     // set the initial maximum value to min possible double
709 jgs 108 return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
710 jgs 102 }
711    
712 jgs 106 double
713     Data::inf() const
714 jgs 102 {
715 jgs 106 //
716     // set the initial minimum value to max possible double
717     return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));
718 jgs 102 }
719    
720     Data
721 jgs 106 Data::maxval() const
722 jgs 102 {
723 jgs 113 //
724     // set the initial maximum value to min possible double
725 jgs 108 return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
726 jgs 102 }
727    
728     Data
729 jgs 106 Data::minval() const
730 jgs 102 {
731 jgs 113 //
732     // set the initial minimum value to max possible double
733 jgs 106 return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));
734 jgs 102 }
735    
736     Data
737 jgs 106 Data::length() const
738 jgs 102 {
739 jgs 106 return dp_algorithm(DataAlgorithmAdapter<Length>(0));
740 jgs 102 }
741    
742     Data
743 jgs 106 Data::trace() const
744 jgs 102 {
745 jgs 106 return dp_algorithm(DataAlgorithmAdapter<Trace>(0));
746 jgs 102 }
747    
748     Data
749 jgs 106 Data::transpose(int axis) const
750 jgs 102 {
751 jgs 106 // not implemented
752     throw DataException("Error - Data::transpose not implemented yet.");
753     return Data();
754 jgs 102 }
755    
756 jgs 104 void
757     Data::saveDX(std::string fileName) const
758     {
759     getDomain().saveDX(fileName,*this);
760     return;
761     }
762    
763 jgs 110 void
764     Data::saveVTK(std::string fileName) const
765     {
766     getDomain().saveVTK(fileName,*this);
767     return;
768     }
769    
770 jgs 102 Data&
771     Data::operator+=(const Data& right)
772     {
773 jgs 94 binaryOp(right,plus<double>());
774     return (*this);
775     }
776    
777 jgs 102 Data&
778     Data::operator+=(const boost::python::object& right)
779 jgs 94 {
780     binaryOp(right,plus<double>());
781     return (*this);
782     }
783    
784 jgs 102 Data&
785     Data::operator-=(const Data& right)
786 jgs 94 {
787     binaryOp(right,minus<double>());
788     return (*this);
789     }
790    
791 jgs 102 Data&
792     Data::operator-=(const boost::python::object& right)
793 jgs 94 {
794     binaryOp(right,minus<double>());
795     return (*this);
796     }
797    
798 jgs 102 Data&
799     Data::operator*=(const Data& right)
800 jgs 94 {
801     binaryOp(right,multiplies<double>());
802     return (*this);
803     }
804    
805 jgs 102 Data&
806     Data::operator*=(const boost::python::object& right)
807 jgs 94 {
808     binaryOp(right,multiplies<double>());
809     return (*this);
810     }
811    
812 jgs 102 Data&
813     Data::operator/=(const Data& right)
814 jgs 94 {
815     binaryOp(right,divides<double>());
816     return (*this);
817     }
818    
819 jgs 102 Data&
820     Data::operator/=(const boost::python::object& right)
821 jgs 94 {
822     binaryOp(right,divides<double>());
823     return (*this);
824     }
825    
826 jgs 102 Data
827     Data::powO(const boost::python::object& right) const
828 jgs 94 {
829     Data result;
830     result.copy(*this);
831     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
832     return result;
833     }
834    
835 jgs 102 Data
836     Data::powD(const Data& right) const
837 jgs 94 {
838     Data result;
839     result.copy(*this);
840     result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
841     return result;
842     }
843    
844     //
845     // NOTE: It is essential to specify the namepsace this operator belongs to
846 jgs 102 Data
847     escript::operator+(const Data& left, const Data& right)
848 jgs 94 {
849     Data result;
850     //
851     // perform a deep copy
852     result.copy(left);
853     result+=right;
854     return result;
855     }
856    
857     //
858     // NOTE: It is essential to specify the namepsace this operator belongs to
859 jgs 102 Data
860     escript::operator-(const Data& left, const Data& right)
861 jgs 94 {
862     Data result;
863     //
864     // perform a deep copy
865     result.copy(left);
866     result-=right;
867     return result;
868     }
869    
870     //
871     // NOTE: It is essential to specify the namepsace this operator belongs to
872 jgs 102 Data
873     escript::operator*(const Data& left, const Data& right)
874 jgs 94 {
875     Data result;
876     //
877     // perform a deep copy
878     result.copy(left);
879     result*=right;
880     return result;
881     }
882    
883     //
884     // NOTE: It is essential to specify the namepsace this operator belongs to
885 jgs 102 Data
886     escript::operator/(const Data& left, const Data& right)
887 jgs 94 {
888     Data result;
889     //
890     // perform a deep copy
891     result.copy(left);
892     result/=right;
893     return result;
894     }
895    
896     //
897     // NOTE: It is essential to specify the namepsace this operator belongs to
898 jgs 102 Data
899     escript::operator+(const Data& left, const boost::python::object& right)
900 jgs 94 {
901     //
902     // Convert to DataArray format if possible
903     DataArray temp(right);
904     Data result;
905     //
906     // perform a deep copy
907     result.copy(left);
908     result+=right;
909     return result;
910     }
911    
912     //
913     // NOTE: It is essential to specify the namepsace this operator belongs to
914 jgs 102 Data
915     escript::operator-(const Data& left, const boost::python::object& right)
916 jgs 94 {
917     //
918     // Convert to DataArray format if possible
919     DataArray temp(right);
920     Data result;
921     //
922     // perform a deep copy
923     result.copy(left);
924     result-=right;
925     return result;
926     }
927    
928     //
929     // NOTE: It is essential to specify the namepsace this operator belongs to
930 jgs 102 Data
931     escript::operator*(const Data& left, const boost::python::object& right)
932 jgs 94 {
933     //
934     // Convert to DataArray format if possible
935     DataArray temp(right);
936     Data result;
937     //
938     // perform a deep copy
939     result.copy(left);
940     result*=right;
941     return result;
942     }
943    
944     //
945     // NOTE: It is essential to specify the namepsace this operator belongs to
946 jgs 102 Data
947     escript::operator/(const Data& left, const boost::python::object& right)
948 jgs 94 {
949     //
950     // Convert to DataArray format if possible
951     DataArray temp(right);
952     Data result;
953     //
954     // perform a deep copy
955     result.copy(left);
956     result/=right;
957     return result;
958     }
959    
960     //
961     // NOTE: It is essential to specify the namepsace this operator belongs to
962 jgs 102 Data
963     escript::operator+(const boost::python::object& left, const Data& right)
964 jgs 94 {
965     //
966     // Construct the result using the given value and the other parameters
967     // from right
968     Data result(left,right);
969     result+=right;
970     return result;
971     }
972    
973     //
974     // NOTE: It is essential to specify the namepsace this operator belongs to
975 jgs 102 Data
976     escript::operator-(const boost::python::object& left, const Data& right)
977 jgs 94 {
978     //
979     // Construct the result using the given value and the other parameters
980     // from right
981     Data result(left,right);
982     result-=right;
983     return result;
984     }
985    
986     //
987     // NOTE: It is essential to specify the namepsace this operator belongs to
988 jgs 102 Data
989     escript::operator*(const boost::python::object& left, const Data& right)
990 jgs 94 {
991     //
992     // Construct the result using the given value and the other parameters
993     // from right
994     Data result(left,right);
995     result*=right;
996     return result;
997     }
998    
999     //
1000     // NOTE: It is essential to specify the namepsace this operator belongs to
1001 jgs 102 Data
1002     escript::operator/(const boost::python::object& left, const Data& right)
1003 jgs 94 {
1004     //
1005     // Construct the result using the given value and the other parameters
1006     // from right
1007     Data result(left,right);
1008     result/=right;
1009     return result;
1010     }
1011    
1012     //
1013     // NOTE: It is essential to specify the namepsace this operator belongs to
1014 jgs 102 //bool escript::operator==(const Data& left, const Data& right)
1015     //{
1016     // /*
1017     // NB: this operator does very little at this point, and isn't to
1018     // be relied on. Requires further implementation.
1019     // */
1020     //
1021     // bool ret;
1022     //
1023     // if (left.isEmpty()) {
1024     // if(!right.isEmpty()) {
1025     // ret = false;
1026     // } else {
1027     // ret = true;
1028     // }
1029     // }
1030     //
1031     // if (left.isConstant()) {
1032     // if(!right.isConstant()) {
1033     // ret = false;
1034     // } else {
1035     // ret = true;
1036     // }
1037     // }
1038     //
1039     // if (left.isTagged()) {
1040     // if(!right.isTagged()) {
1041     // ret = false;
1042     // } else {
1043     // ret = true;
1044     // }
1045     // }
1046     //
1047     // if (left.isExpanded()) {
1048     // if(!right.isExpanded()) {
1049     // ret = false;
1050     // } else {
1051     // ret = true;
1052     // }
1053     // }
1054     //
1055     // return ret;
1056     //}
1057    
1058     Data
1059     Data::getItem(const boost::python::object& key) const
1060 jgs 94 {
1061 jgs 102 const DataArrayView& view=getPointDataView();
1062 jgs 94
1063 jgs 102 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1064 jgs 94
1065 jgs 102 if (slice_region.size()!=view.getRank()) {
1066     throw DataException("Error - slice size does not match Data rank.");
1067 jgs 94 }
1068    
1069 jgs 102 return getSlice(slice_region);
1070 jgs 94 }
1071    
1072     Data
1073 jgs 102 Data::getSlice(const DataArrayView::RegionType& region) const
1074 jgs 94 {
1075 jgs 102 return Data(*this,region);
1076 jgs 94 }
1077    
1078     void
1079 jgs 102 Data::setItemO(const boost::python::object& key,
1080     const boost::python::object& value)
1081 jgs 94 {
1082 jgs 102 Data tempData(value,getFunctionSpace());
1083     setItemD(key,tempData);
1084     }
1085    
1086     void
1087     Data::setItemD(const boost::python::object& key,
1088     const Data& value)
1089     {
1090 jgs 94 const DataArrayView& view=getPointDataView();
1091     DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1092     if (slice_region.size()!=view.getRank()) {
1093     throw DataException("Error - slice size does not match Data rank.");
1094     }
1095 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
1096     setSlice(Data(value,getFunctionSpace()),slice_region);
1097     } else {
1098     setSlice(value,slice_region);
1099     }
1100 jgs 94 }
1101    
1102     void
1103 jgs 102 Data::setSlice(const Data& value,
1104     const DataArrayView::RegionType& region)
1105 jgs 94 {
1106 jgs 102 Data tempValue(value);
1107     typeMatchLeft(tempValue);
1108     typeMatchRight(tempValue);
1109     m_data->setSlice(tempValue.m_data.get(),region);
1110     }
1111    
1112     void
1113     Data::typeMatchLeft(Data& right) const
1114     {
1115     if (isExpanded()){
1116     right.expand();
1117     } else if (isTagged()) {
1118     if (right.isConstant()) {
1119     right.tag();
1120     }
1121     }
1122     }
1123    
1124     void
1125     Data::typeMatchRight(const Data& right)
1126     {
1127 jgs 94 if (isTagged()) {
1128     if (right.isExpanded()) {
1129     expand();
1130     }
1131     } else if (isConstant()) {
1132     if (right.isExpanded()) {
1133     expand();
1134     } else if (right.isTagged()) {
1135     tag();
1136     }
1137     }
1138     }
1139    
1140     void
1141     Data::setTaggedValue(int tagKey,
1142     const boost::python::object& value)
1143     {
1144     //
1145     // Ensure underlying data object is of type DataTagged
1146     tag();
1147    
1148     if (!isTagged()) {
1149     throw DataException("Error - DataTagged conversion failed!!");
1150     }
1151    
1152     //
1153     // Construct DataArray from boost::python::object input value
1154     DataArray valueDataArray(value);
1155    
1156     //
1157     // Call DataAbstract::setTaggedValue
1158     m_data->setTaggedValue(tagKey,valueDataArray.getView());
1159     }
1160    
1161 jgs 110 void
1162     Data::setRefValue(int ref,
1163     const boost::python::numeric::array& value)
1164     {
1165     //
1166     // Construct DataArray from boost::python::object input value
1167     DataArray valueDataArray(value);
1168    
1169     //
1170     // Call DataAbstract::setRefValue
1171     m_data->setRefValue(ref,valueDataArray);
1172     }
1173    
1174     void
1175     Data::getRefValue(int ref,
1176     boost::python::numeric::array& value)
1177     {
1178     //
1179     // Construct DataArray for boost::python::object return value
1180     DataArray valueDataArray(value);
1181    
1182     //
1183     // Load DataArray with values from data-points specified by ref
1184     m_data->getRefValue(ref,valueDataArray);
1185    
1186     //
1187     // Load values from valueDataArray into return numarray
1188    
1189     // extract the shape of the numarray
1190     int rank = value.getrank();
1191     DataArrayView::ShapeType shape;
1192     for (int i=0; i < rank; i++) {
1193     shape.push_back(extract<int>(value.getshape()[i]));
1194     }
1195    
1196     // and load the numarray with the data from the DataArray
1197     DataArrayView valueView = valueDataArray.getView();
1198    
1199     if (rank==0) {
1200     throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1201     }
1202     if (rank==1) {
1203     for (int i=0; i < shape[0]; i++) {
1204     value[i] = valueView(i);
1205     }
1206     }
1207     if (rank==2) {
1208     throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1209     }
1210     if (rank==3) {
1211     throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1212     }
1213     if (rank==4) {
1214     throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1215     }
1216    
1217     }
1218    
1219 jgs 94 /*
1220     Note: this version removed for now. Not needed, and breaks escript.cpp
1221     void
1222     Data::setTaggedValue(int tagKey,
1223     const DataArrayView& value)
1224     {
1225     //
1226     // Ensure underlying data object is of type DataTagged
1227     tag();
1228    
1229     if (!isTagged()) {
1230     throw DataException("Error - DataTagged conversion failed!!");
1231     }
1232    
1233     //
1234     // Call DataAbstract::setTaggedValue
1235     m_data->setTaggedValue(tagKey,value);
1236     }
1237     */
1238    
1239 jgs 119 void
1240     Data::archiveData(const std::string fileName)
1241     {
1242     cout << "Archiving Data object to: " << fileName << endl;
1243    
1244     //
1245     // Determine type of this Data object
1246     int dataType = -1;
1247    
1248     if (isEmpty()) {
1249     dataType = 0;
1250     cout << "\tdataType: DataEmpty" << endl;
1251     }
1252     if (isConstant()) {
1253     dataType = 1;
1254     cout << "\tdataType: DataConstant" << endl;
1255     }
1256     if (isTagged()) {
1257     dataType = 2;
1258     cout << "\tdataType: DataTagged" << endl;
1259     }
1260     if (isExpanded()) {
1261     dataType = 3;
1262     cout << "\tdataType: DataExpanded" << endl;
1263     }
1264     if (dataType == -1) {
1265     throw DataException("archiveData Error: undefined dataType");
1266     }
1267    
1268     //
1269     // Collect data items common to all Data types
1270     int noSamples = getNumSamples();
1271     int noDPPSample = getNumDataPointsPerSample();
1272     int functionSpaceType = getFunctionSpace().getTypeCode();
1273     int dataPointRank = getDataPointRank();
1274     int dataPointSize = getDataPointSize();
1275     int dataLength = getLength();
1276     DataArrayView::ShapeType dataPointShape = getDataPointShape();
1277     int referenceNumbers[noSamples];
1278     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1279     referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1280     }
1281     int tagNumbers[noSamples];
1282     if (isTagged()) {
1283     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1284     tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1285     }
1286     }
1287    
1288     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1289     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1290     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1291    
1292     //
1293     // Flatten Shape to an array of integers suitable for writing to file
1294     int flatShape[4] = {0,0,0,0};
1295     cout << "\tshape: < ";
1296     for (int dim=0; dim<dataPointRank; dim++) {
1297     flatShape[dim] = dataPointShape[dim];
1298     cout << dataPointShape[dim] << " ";
1299     }
1300     cout << ">" << endl;
1301    
1302     //
1303     // Write common data items to archive file
1304     ofstream archiveFile;
1305     archiveFile.open(fileName.data(), ios::out);
1306    
1307     if (!archiveFile.good()) {
1308     throw DataException("archiveData Error: problem opening archive file");
1309     }
1310    
1311     archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1312     archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1313     archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1314     archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1315     archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1316     archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1317     archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1318     for (int dim = 0; dim < 4; dim++) {
1319     archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1320     }
1321     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1322     archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1323     }
1324     if (isTagged()) {
1325     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1326     archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1327     }
1328     }
1329    
1330     if (!archiveFile.good()) {
1331     throw DataException("archiveData Error: problem writing to archive file");
1332     }
1333    
1334     archiveFile.close();
1335    
1336     if (!archiveFile.good()) {
1337     throw DataException("archiveData Error: problem closing archive file");
1338     }
1339    
1340     //
1341     // Collect and archive underlying data values for each Data type
1342     switch (dataType) {
1343     case 0:
1344     // DataEmpty
1345     break;
1346     case 1:
1347     // DataConstant
1348     break;
1349     case 2:
1350     // DataTagged
1351     break;
1352     case 3:
1353     // DataExpanded
1354     break;
1355     }
1356    
1357     }
1358    
1359     void
1360     Data::extractData(const std::string fileName,
1361     const FunctionSpace& fspace)
1362     {
1363     //
1364     // Can only extract Data to an object which is initially DataEmpty
1365     if (!isEmpty()) {
1366     throw DataException("extractData Error: can only extract to DataEmpty object");
1367     }
1368    
1369     cout << "Extracting Data object from: " << fileName << endl;
1370    
1371     int dataType;
1372     int noSamples;
1373     int noDPPSample;
1374     int functionSpaceType;
1375     int dataPointRank;
1376     int dataPointSize;
1377     int dataLength;
1378     DataArrayView::ShapeType dataPointShape;
1379     int flatShape[4];
1380    
1381     //
1382     // Open the archive file and read common data items
1383     ifstream archiveFile;
1384     archiveFile.open(fileName.data(), ios::in);
1385    
1386     if (!archiveFile.good()) {
1387     throw DataException("extractData Error: problem opening archive file");
1388     }
1389    
1390     archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1391     archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1392     archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1393     archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1394     archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1395     archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1396     archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1397     for (int dim = 0; dim < 4; dim++) {
1398     archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1399     if (flatShape[dim]>0) {
1400     dataPointShape.push_back(flatShape[dim]);
1401     }
1402     }
1403     int referenceNumbers[noSamples];
1404     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1405     archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1406     }
1407     int tagNumbers[noSamples];
1408     if (dataType==2) {
1409     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1410     archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1411     }
1412     }
1413    
1414     if (!archiveFile.good()) {
1415     throw DataException("extractData Error: problem reading from archive file");
1416     }
1417    
1418     archiveFile.close();
1419    
1420     if (!archiveFile.good()) {
1421     throw DataException("extractData Error: problem closing archive file");
1422     }
1423    
1424     switch (dataType) {
1425     case 0:
1426     cout << "\tdataType: DataEmpty" << endl;
1427     break;
1428     case 1:
1429     cout << "\tdataType: DataConstant" << endl;
1430     break;
1431     case 2:
1432     cout << "\tdataType: DataTagged" << endl;
1433     break;
1434     case 3:
1435     cout << "\tdataType: DataExpanded" << endl;
1436     break;
1437     default:
1438     throw DataException("extractData Error: undefined dataType read from archive file");
1439     break;
1440     }
1441    
1442     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1443     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1444     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1445     cout << "\tshape: < ";
1446     for (int dim = 0; dim < dataPointRank; dim++) {
1447     cout << dataPointShape[dim] << " ";
1448     }
1449     cout << ">" << endl;
1450    
1451     //
1452     // Verify that supplied FunctionSpace object is compatible with this Data object.
1453     if ( (fspace.getTypeCode()!=functionSpaceType) ||
1454     (fspace.getNumSamples()!=noSamples) ||
1455     (fspace.getNumDPPSample()!=noDPPSample)
1456     ) {
1457     throw DataException("extractData Error: incompatible FunctionSpace");
1458     }
1459     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1460     if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
1461     throw DataException("extractData Error: incompatible FunctionSpace");
1462     }
1463     }
1464     if (dataType==2) {
1465     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1466     if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
1467     throw DataException("extractData Error: incompatible FunctionSpace");
1468     }
1469     }
1470     }
1471    
1472     //
1473     // Construct a DataVector to hold underlying data values
1474     DataVector dataVec(dataLength);
1475    
1476     //
1477     // Load this DataVector with the appropriate values
1478     switch (dataType) {
1479     case 0:
1480     // DataEmpty
1481     break;
1482     case 1:
1483     // DataConstant
1484     break;
1485     case 2:
1486     // DataTagged
1487     break;
1488     case 3:
1489     // DataExpanded
1490     break;
1491     }
1492    
1493     //
1494     // Construct an appropriate Data object
1495     DataAbstract* tempData;
1496     switch (dataType) {
1497     case 0:
1498     // DataEmpty
1499     tempData=new DataEmpty();
1500     break;
1501     case 1:
1502     // DataConstant
1503     tempData=new DataConstant(fspace,dataPointShape,dataVec);
1504     break;
1505     case 2:
1506     // DataTagged
1507     tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
1508     break;
1509     case 3:
1510     // DataExpanded
1511     tempData=new DataExpanded(fspace,dataPointShape,dataVec);
1512     break;
1513     }
1514     shared_ptr<DataAbstract> temp_data(tempData);
1515     m_data=temp_data;
1516    
1517     }
1518    
1519 jgs 94 ostream& escript::operator<<(ostream& o, const Data& data)
1520     {
1521     o << data.toString();
1522     return o;
1523     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26