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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26