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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 151 - (show annotations)
Thu Sep 22 01:55:00 2005 UTC (14 years ago) by jgs
Original Path: trunk/esys2/escript/src/Data/Data.cpp
File size: 52056 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-22

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26