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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 480 - (show annotations)
Wed Feb 1 05:15:12 2006 UTC (13 years, 7 months ago) by jgs
File size: 51939 byte(s)
rationalise #includes and forward declarations

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26