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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26