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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (show annotations)
Wed Nov 9 02:02:19 2005 UTC (14 years, 1 month ago) by jgs
Original Path: trunk/escript/src/Data/Data.cpp
File size: 52231 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26