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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 117 - (show annotations)
Fri Apr 1 05:48:57 2005 UTC (14 years, 7 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/Data.cpp
File size: 28877 byte(s)
*** empty log message ***

1 // $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 #include <algorithm>
22 #include <vector>
23 #include <exception>
24 #include <functional>
25 #include <math.h>
26
27 #include <boost/python/str.hpp>
28 #include <boost/python/extract.hpp>
29 #include <boost/python/long.hpp>
30
31 #include "escript/Data/DataException.h"
32
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 #include "escript/Data/UnaryFuncs.h"
42
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 shared_ptr<DataAbstract> temp_data(temp);
54 m_data=temp_data;
55 }
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 Data::Data(const Data& inData)
81 {
82 m_data=inData.m_data;
83 }
84
85 Data::Data(const Data& inData,
86 const DataArrayView::RegionType& region)
87 {
88 //
89 // 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 }
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 shared_ptr<DataAbstract> temp_data(temp);
124 m_data=temp_data;
125 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 shared_ptr<DataAbstract> temp_data(newData);
220 m_data=temp_data;
221 return;
222 }
223 }
224 {
225 DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());
226 if (temp!=0) {
227 //
228 // Construct a DataTagged copy
229 DataAbstract* newData=new DataTagged(*temp);
230 shared_ptr<DataAbstract> temp_data(newData);
231 m_data=temp_data;
232 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 shared_ptr<DataAbstract> temp_data(newData);
242 m_data=temp_data;
243 return;
244 }
245 }
246 {
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 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 shared_ptr<DataAbstract> temp_data(temp);
312 m_data=temp_data;
313 } else if (isTagged()) {
314 DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
315 DataAbstract* temp=new DataExpanded(*tempDataTag);
316 shared_ptr<DataAbstract> temp_data(temp);
317 m_data=temp_data;
318 } 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 shared_ptr<DataAbstract> temp_data(temp);
335 m_data=temp_data;
336 } 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 void
348 Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)
349 {
350 m_data->reshapeDataPoint(shape);
351 }
352
353 Data
354 Data::wherePositive() const
355 {
356 return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
357 }
358
359 Data
360 Data::whereNegative() const
361 {
362 return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
363 }
364
365 Data
366 Data::whereNonNegative() const
367 {
368 return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
369 }
370
371 Data
372 Data::whereNonPositive() const
373 {
374 return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
375 }
376
377 Data
378 Data::whereZero() const
379 {
380 return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));
381 }
382
383 Data
384 Data::whereNonZero() const
385 {
386 return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));
387 }
388
389 Data
390 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 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 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 bp_array.resize(1);
573 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 Data
653 Data::sign() const
654 {
655 return escript::unaryOp(*this,escript::fsign);
656 }
657
658 Data
659 Data::abs() const
660 {
661 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
662 }
663
664 Data
665 Data::neg() const
666 {
667 return escript::unaryOp(*this,negate<double>());
668 }
669
670 Data
671 Data::pos() const
672 {
673 return (*this);
674 }
675
676 Data
677 Data::exp() const
678 {
679 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
680 }
681
682 Data
683 Data::sqrt() const
684 {
685 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
686 }
687
688 double
689 Data::Lsup() const
690 {
691 //
692 // set the initial absolute maximum value to zero
693 return algorithm(DataAlgorithmAdapter<AbsMax>(0));
694 }
695
696 double
697 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 Data::sup() const
706 {
707 //
708 // set the initial maximum value to min possible double
709 return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
710 }
711
712 double
713 Data::inf() const
714 {
715 //
716 // set the initial minimum value to max possible double
717 return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));
718 }
719
720 Data
721 Data::maxval() const
722 {
723 //
724 // set the initial maximum value to min possible double
725 return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
726 }
727
728 Data
729 Data::minval() const
730 {
731 //
732 // set the initial minimum value to max possible double
733 return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));
734 }
735
736 Data
737 Data::length() const
738 {
739 return dp_algorithm(DataAlgorithmAdapter<Length>(0));
740 }
741
742 Data
743 Data::trace() const
744 {
745 return dp_algorithm(DataAlgorithmAdapter<Trace>(0));
746 }
747
748 Data
749 Data::transpose(int axis) const
750 {
751 // not implemented
752 throw DataException("Error - Data::transpose not implemented yet.");
753 return Data();
754 }
755
756 void
757 Data::saveDX(std::string fileName) const
758 {
759 getDomain().saveDX(fileName,*this);
760 return;
761 }
762
763 void
764 Data::saveVTK(std::string fileName) const
765 {
766 getDomain().saveVTK(fileName,*this);
767 return;
768 }
769
770 Data&
771 Data::operator+=(const Data& right)
772 {
773 binaryOp(right,plus<double>());
774 return (*this);
775 }
776
777 Data&
778 Data::operator+=(const boost::python::object& right)
779 {
780 binaryOp(right,plus<double>());
781 return (*this);
782 }
783
784 Data&
785 Data::operator-=(const Data& right)
786 {
787 binaryOp(right,minus<double>());
788 return (*this);
789 }
790
791 Data&
792 Data::operator-=(const boost::python::object& right)
793 {
794 binaryOp(right,minus<double>());
795 return (*this);
796 }
797
798 Data&
799 Data::operator*=(const Data& right)
800 {
801 binaryOp(right,multiplies<double>());
802 return (*this);
803 }
804
805 Data&
806 Data::operator*=(const boost::python::object& right)
807 {
808 binaryOp(right,multiplies<double>());
809 return (*this);
810 }
811
812 Data&
813 Data::operator/=(const Data& right)
814 {
815 binaryOp(right,divides<double>());
816 return (*this);
817 }
818
819 Data&
820 Data::operator/=(const boost::python::object& right)
821 {
822 binaryOp(right,divides<double>());
823 return (*this);
824 }
825
826 Data
827 Data::powO(const boost::python::object& right) const
828 {
829 Data result;
830 result.copy(*this);
831 result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
832 return result;
833 }
834
835 Data
836 Data::powD(const Data& right) const
837 {
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 Data
847 escript::operator+(const Data& left, const Data& right)
848 {
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 Data
860 escript::operator-(const Data& left, const Data& right)
861 {
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 Data
873 escript::operator*(const Data& left, const Data& right)
874 {
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 Data
886 escript::operator/(const Data& left, const Data& right)
887 {
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 Data
899 escript::operator+(const Data& left, const boost::python::object& right)
900 {
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 Data
915 escript::operator-(const Data& left, const boost::python::object& right)
916 {
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 Data
931 escript::operator*(const Data& left, const boost::python::object& right)
932 {
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 Data
947 escript::operator/(const Data& left, const boost::python::object& right)
948 {
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 Data
963 escript::operator+(const boost::python::object& left, const Data& right)
964 {
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 Data
976 escript::operator-(const boost::python::object& left, const Data& right)
977 {
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 Data
989 escript::operator*(const boost::python::object& left, const Data& right)
990 {
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 Data
1002 escript::operator/(const boost::python::object& left, const Data& right)
1003 {
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 //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 {
1061 const DataArrayView& view=getPointDataView();
1062
1063 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1064
1065 if (slice_region.size()!=view.getRank()) {
1066 throw DataException("Error - slice size does not match Data rank.");
1067 }
1068
1069 return getSlice(slice_region);
1070 }
1071
1072 Data
1073 Data::getSlice(const DataArrayView::RegionType& region) const
1074 {
1075 return Data(*this,region);
1076 }
1077
1078 void
1079 Data::setItemO(const boost::python::object& key,
1080 const boost::python::object& value)
1081 {
1082 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 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 if (getFunctionSpace()!=value.getFunctionSpace()) {
1096 setSlice(Data(value,getFunctionSpace()),slice_region);
1097 } else {
1098 setSlice(value,slice_region);
1099 }
1100 }
1101
1102 void
1103 Data::setSlice(const Data& value,
1104 const DataArrayView::RegionType& region)
1105 {
1106 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 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 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 /*
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 ostream& escript::operator<<(ostream& o, const Data& data)
1240 {
1241 o << data.toString();
1242 return o;
1243 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26