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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 436 - (show annotations)
Thu Jan 19 22:36:36 2006 UTC (13 years, 10 months ago) by gross
Original Path: trunk/escript/src/Data/Data.cpp
File size: 52188 byte(s)
bug fixed in integrate method. don;t use a[i,j]=.. to set a value in the numeric::array a but 
a[make_tuple(i,j)].


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26