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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 615 - (show annotations)
Wed Mar 22 02:12:00 2006 UTC (13 years, 6 months ago) by elspeth
File size: 53199 byte(s)
More copyright information.

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(double tol) const
405 {
406 #if defined DOPROF
407 profData->where++;
408 #endif
409 return escript::unaryOp(*this,bind2nd(greater<double>(),0.0+tol));
410 }
411
412 Data
413 Data::whereNegative(double tol) const
414 {
415 #if defined DOPROF
416 profData->where++;
417 #endif
418 return escript::unaryOp(*this,bind2nd(less<double>(),0.0-tol));
419 }
420
421 Data
422 Data::whereNonNegative(double tol) const
423 {
424 #if defined DOPROF
425 profData->where++;
426 #endif
427 return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0+tol));
428 }
429
430 Data
431 Data::whereNonPositive(double tol) const
432 {
433 #if defined DOPROF
434 profData->where++;
435 #endif
436 return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0-tol));
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::powO(const boost::python::object& right) const
1393 {
1394 #if defined DOPROF
1395 profData->binary++;
1396 #endif
1397 Data result;
1398 result.copy(*this);
1399 result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1400 return result;
1401 }
1402
1403 Data
1404 Data::powD(const Data& right) const
1405 {
1406 #if defined DOPROF
1407 profData->binary++;
1408 #endif
1409 Data result;
1410 result.copy(*this);
1411 result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1412 return result;
1413 }
1414
1415 //
1416 // NOTE: It is essential to specify the namespace this operator belongs to
1417 Data
1418 escript::operator+(const Data& left, const Data& right)
1419 {
1420 Data result;
1421 //
1422 // perform a deep copy
1423 result.copy(left);
1424 result+=right;
1425 return result;
1426 }
1427
1428 //
1429 // NOTE: It is essential to specify the namespace this operator belongs to
1430 Data
1431 escript::operator-(const Data& left, const Data& right)
1432 {
1433 Data result;
1434 //
1435 // perform a deep copy
1436 result.copy(left);
1437 result-=right;
1438 return result;
1439 }
1440
1441 //
1442 // NOTE: It is essential to specify the namespace this operator belongs to
1443 Data
1444 escript::operator*(const Data& left, const Data& right)
1445 {
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 Data& right)
1458 {
1459 Data result;
1460 //
1461 // perform a deep copy
1462 result.copy(left);
1463 result/=right;
1464 return result;
1465 }
1466
1467 //
1468 // NOTE: It is essential to specify the namespace this operator belongs to
1469 Data
1470 escript::operator+(const Data& left, const boost::python::object& right)
1471 {
1472 //
1473 // Convert to DataArray format if possible
1474 DataArray temp(right);
1475 Data result;
1476 //
1477 // perform a deep copy
1478 result.copy(left);
1479 result+=right;
1480 return result;
1481 }
1482
1483 //
1484 // NOTE: It is essential to specify the namespace this operator belongs to
1485 Data
1486 escript::operator-(const Data& left, const boost::python::object& right)
1487 {
1488 //
1489 // Convert to DataArray format if possible
1490 DataArray temp(right);
1491 Data result;
1492 //
1493 // perform a deep copy
1494 result.copy(left);
1495 result-=right;
1496 return result;
1497 }
1498
1499 //
1500 // NOTE: It is essential to specify the namespace this operator belongs to
1501 Data
1502 escript::operator*(const Data& left, const boost::python::object& right)
1503 {
1504 //
1505 // Convert to DataArray format if possible
1506 DataArray temp(right);
1507 Data result;
1508 //
1509 // perform a deep copy
1510 result.copy(left);
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 Data& left, const boost::python::object& right)
1519 {
1520 //
1521 // Convert to DataArray format if possible
1522 DataArray temp(right);
1523 Data result;
1524 //
1525 // perform a deep copy
1526 result.copy(left);
1527 result/=right;
1528 return result;
1529 }
1530
1531 //
1532 // NOTE: It is essential to specify the namespace this operator belongs to
1533 Data
1534 escript::operator+(const boost::python::object& left, const Data& right)
1535 {
1536 //
1537 // Construct the result using the given value and the other parameters
1538 // from right
1539 Data result(left,right);
1540 result+=right;
1541 return result;
1542 }
1543
1544 //
1545 // NOTE: It is essential to specify the namespace this operator belongs to
1546 Data
1547 escript::operator-(const boost::python::object& left, const Data& right)
1548 {
1549 //
1550 // Construct the result using the given value and the other parameters
1551 // from right
1552 Data result(left,right);
1553 result-=right;
1554 return result;
1555 }
1556
1557 //
1558 // NOTE: It is essential to specify the namespace this operator belongs to
1559 Data
1560 escript::operator*(const boost::python::object& left, const Data& right)
1561 {
1562 //
1563 // Construct the result using the given value and the other parameters
1564 // from right
1565 Data result(left,right);
1566 result*=right;
1567 return result;
1568 }
1569
1570 //
1571 // NOTE: It is essential to specify the namespace this operator belongs to
1572 Data
1573 escript::operator/(const boost::python::object& left, const Data& right)
1574 {
1575 //
1576 // Construct the result using the given value and the other parameters
1577 // from right
1578 Data result(left,right);
1579 result/=right;
1580 return result;
1581 }
1582
1583 //
1584 //bool escript::operator==(const Data& left, const Data& right)
1585 //{
1586 // /*
1587 // NB: this operator does very little at this point, and isn't to
1588 // be relied on. Requires further implementation.
1589 // */
1590 //
1591 // bool ret;
1592 //
1593 // if (left.isEmpty()) {
1594 // if(!right.isEmpty()) {
1595 // ret = false;
1596 // } else {
1597 // ret = true;
1598 // }
1599 // }
1600 //
1601 // if (left.isConstant()) {
1602 // if(!right.isConstant()) {
1603 // ret = false;
1604 // } else {
1605 // ret = true;
1606 // }
1607 // }
1608 //
1609 // if (left.isTagged()) {
1610 // if(!right.isTagged()) {
1611 // ret = false;
1612 // } else {
1613 // ret = true;
1614 // }
1615 // }
1616 //
1617 // if (left.isExpanded()) {
1618 // if(!right.isExpanded()) {
1619 // ret = false;
1620 // } else {
1621 // ret = true;
1622 // }
1623 // }
1624 //
1625 // return ret;
1626 //}
1627
1628 Data
1629 Data::getItem(const boost::python::object& key) const
1630 {
1631 const DataArrayView& view=getPointDataView();
1632
1633 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1634
1635 if (slice_region.size()!=view.getRank()) {
1636 throw DataException("Error - slice size does not match Data rank.");
1637 }
1638
1639 return getSlice(slice_region);
1640 }
1641
1642 Data
1643 Data::getSlice(const DataArrayView::RegionType& region) const
1644 {
1645 #if defined DOPROF
1646 profData->slicing++;
1647 #endif
1648 return Data(*this,region);
1649 }
1650
1651 void
1652 Data::setItemO(const boost::python::object& key,
1653 const boost::python::object& value)
1654 {
1655 Data tempData(value,getFunctionSpace());
1656 setItemD(key,tempData);
1657 }
1658
1659 void
1660 Data::setItemD(const boost::python::object& key,
1661 const Data& value)
1662 {
1663 const DataArrayView& view=getPointDataView();
1664
1665 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1666 if (slice_region.size()!=view.getRank()) {
1667 throw DataException("Error - slice size does not match Data rank.");
1668 }
1669 if (getFunctionSpace()!=value.getFunctionSpace()) {
1670 setSlice(Data(value,getFunctionSpace()),slice_region);
1671 } else {
1672 setSlice(value,slice_region);
1673 }
1674 }
1675
1676 void
1677 Data::setSlice(const Data& value,
1678 const DataArrayView::RegionType& region)
1679 {
1680 #if defined DOPROF
1681 profData->slicing++;
1682 #endif
1683 Data tempValue(value);
1684 typeMatchLeft(tempValue);
1685 typeMatchRight(tempValue);
1686 m_data->setSlice(tempValue.m_data.get(),region);
1687 }
1688
1689 void
1690 Data::typeMatchLeft(Data& right) const
1691 {
1692 if (isExpanded()){
1693 right.expand();
1694 } else if (isTagged()) {
1695 if (right.isConstant()) {
1696 right.tag();
1697 }
1698 }
1699 }
1700
1701 void
1702 Data::typeMatchRight(const Data& right)
1703 {
1704 if (isTagged()) {
1705 if (right.isExpanded()) {
1706 expand();
1707 }
1708 } else if (isConstant()) {
1709 if (right.isExpanded()) {
1710 expand();
1711 } else if (right.isTagged()) {
1712 tag();
1713 }
1714 }
1715 }
1716
1717 void
1718 Data::setTaggedValue(int tagKey,
1719 const boost::python::object& value)
1720 {
1721 //
1722 // Ensure underlying data object is of type DataTagged
1723 tag();
1724
1725 if (!isTagged()) {
1726 throw DataException("Error - DataTagged conversion failed!!");
1727 }
1728
1729 //
1730 // Construct DataArray from boost::python::object input value
1731 DataArray valueDataArray(value);
1732
1733 //
1734 // Call DataAbstract::setTaggedValue
1735 m_data->setTaggedValue(tagKey,valueDataArray.getView());
1736 }
1737
1738 void
1739 Data::setTaggedValueFromCPP(int tagKey,
1740 const DataArrayView& value)
1741 {
1742 //
1743 // Ensure underlying data object is of type DataTagged
1744 tag();
1745
1746 if (!isTagged()) {
1747 throw DataException("Error - DataTagged conversion failed!!");
1748 }
1749
1750 //
1751 // Call DataAbstract::setTaggedValue
1752 m_data->setTaggedValue(tagKey,value);
1753 }
1754
1755 int
1756 Data::getTagNumber(int dpno)
1757 {
1758 return m_data->getTagNumber(dpno);
1759 }
1760
1761 void
1762 Data::setRefValue(int ref,
1763 const boost::python::numeric::array& value)
1764 {
1765 //
1766 // Construct DataArray from boost::python::object input value
1767 DataArray valueDataArray(value);
1768
1769 //
1770 // Call DataAbstract::setRefValue
1771 m_data->setRefValue(ref,valueDataArray);
1772 }
1773
1774 void
1775 Data::getRefValue(int ref,
1776 boost::python::numeric::array& value)
1777 {
1778 //
1779 // Construct DataArray for boost::python::object return value
1780 DataArray valueDataArray(value);
1781
1782 //
1783 // Load DataArray with values from data-points specified by ref
1784 m_data->getRefValue(ref,valueDataArray);
1785
1786 //
1787 // Load values from valueDataArray into return numarray
1788
1789 // extract the shape of the numarray
1790 int rank = value.getrank();
1791 DataArrayView::ShapeType shape;
1792 for (int i=0; i < rank; i++) {
1793 shape.push_back(extract<int>(value.getshape()[i]));
1794 }
1795
1796 // and load the numarray with the data from the DataArray
1797 DataArrayView valueView = valueDataArray.getView();
1798
1799 if (rank==0) {
1800 boost::python::numeric::array temp_numArray(valueView());
1801 value = temp_numArray;
1802 }
1803 if (rank==1) {
1804 for (int i=0; i < shape[0]; i++) {
1805 value[i] = valueView(i);
1806 }
1807 }
1808 if (rank==2) {
1809 for (int i=0; i < shape[0]; i++) {
1810 for (int j=0; j < shape[1]; j++) {
1811 value[i][j] = valueView(i,j);
1812 }
1813 }
1814 }
1815 if (rank==3) {
1816 for (int i=0; i < shape[0]; i++) {
1817 for (int j=0; j < shape[1]; j++) {
1818 for (int k=0; k < shape[2]; k++) {
1819 value[i][j][k] = valueView(i,j,k);
1820 }
1821 }
1822 }
1823 }
1824 if (rank==4) {
1825 for (int i=0; i < shape[0]; i++) {
1826 for (int j=0; j < shape[1]; j++) {
1827 for (int k=0; k < shape[2]; k++) {
1828 for (int l=0; l < shape[3]; l++) {
1829 value[i][j][k][l] = valueView(i,j,k,l);
1830 }
1831 }
1832 }
1833 }
1834 }
1835
1836 }
1837
1838 void
1839 Data::archiveData(const std::string fileName)
1840 {
1841 cout << "Archiving Data object to: " << fileName << endl;
1842
1843 //
1844 // Determine type of this Data object
1845 int dataType = -1;
1846
1847 if (isEmpty()) {
1848 dataType = 0;
1849 cout << "\tdataType: DataEmpty" << endl;
1850 }
1851 if (isConstant()) {
1852 dataType = 1;
1853 cout << "\tdataType: DataConstant" << endl;
1854 }
1855 if (isTagged()) {
1856 dataType = 2;
1857 cout << "\tdataType: DataTagged" << endl;
1858 }
1859 if (isExpanded()) {
1860 dataType = 3;
1861 cout << "\tdataType: DataExpanded" << endl;
1862 }
1863
1864 if (dataType == -1) {
1865 throw DataException("archiveData Error: undefined dataType");
1866 }
1867
1868 //
1869 // Collect data items common to all Data types
1870 int noSamples = getNumSamples();
1871 int noDPPSample = getNumDataPointsPerSample();
1872 int functionSpaceType = getFunctionSpace().getTypeCode();
1873 int dataPointRank = getDataPointRank();
1874 int dataPointSize = getDataPointSize();
1875 int dataLength = getLength();
1876 DataArrayView::ShapeType dataPointShape = getDataPointShape();
1877 int referenceNumbers[noSamples];
1878 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1879 referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1880 }
1881 int tagNumbers[noSamples];
1882 if (isTagged()) {
1883 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1884 tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1885 }
1886 }
1887
1888 cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1889 cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1890 cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1891
1892 //
1893 // Flatten Shape to an array of integers suitable for writing to file
1894 int flatShape[4] = {0,0,0,0};
1895 cout << "\tshape: < ";
1896 for (int dim=0; dim<dataPointRank; dim++) {
1897 flatShape[dim] = dataPointShape[dim];
1898 cout << dataPointShape[dim] << " ";
1899 }
1900 cout << ">" << endl;
1901
1902 //
1903 // Open archive file
1904 ofstream archiveFile;
1905 archiveFile.open(fileName.data(), ios::out);
1906
1907 if (!archiveFile.good()) {
1908 throw DataException("archiveData Error: problem opening archive file");
1909 }
1910
1911 //
1912 // Write common data items to archive file
1913 archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1914 archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1915 archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1916 archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1917 archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1918 archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1919 archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1920 for (int dim = 0; dim < 4; dim++) {
1921 archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1922 }
1923 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1924 archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1925 }
1926 if (isTagged()) {
1927 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1928 archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1929 }
1930 }
1931
1932 if (!archiveFile.good()) {
1933 throw DataException("archiveData Error: problem writing to archive file");
1934 }
1935
1936 //
1937 // Archive underlying data values for each Data type
1938 int noValues;
1939 switch (dataType) {
1940 case 0:
1941 // DataEmpty
1942 noValues = 0;
1943 archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1944 cout << "\tnoValues: " << noValues << endl;
1945 break;
1946 case 1:
1947 // DataConstant
1948 noValues = m_data->getLength();
1949 archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1950 cout << "\tnoValues: " << noValues << endl;
1951 if (m_data->archiveData(archiveFile,noValues)) {
1952 throw DataException("archiveData Error: problem writing data to archive file");
1953 }
1954 break;
1955 case 2:
1956 // DataTagged
1957 noValues = m_data->getLength();
1958 archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1959 cout << "\tnoValues: " << noValues << endl;
1960 if (m_data->archiveData(archiveFile,noValues)) {
1961 throw DataException("archiveData Error: problem writing data to archive file");
1962 }
1963 break;
1964 case 3:
1965 // DataExpanded
1966 noValues = m_data->getLength();
1967 archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1968 cout << "\tnoValues: " << noValues << endl;
1969 if (m_data->archiveData(archiveFile,noValues)) {
1970 throw DataException("archiveData Error: problem writing data to archive file");
1971 }
1972 break;
1973 }
1974
1975 if (!archiveFile.good()) {
1976 throw DataException("archiveData Error: problem writing data to archive file");
1977 }
1978
1979 //
1980 // Close archive file
1981 archiveFile.close();
1982
1983 if (!archiveFile.good()) {
1984 throw DataException("archiveData Error: problem closing archive file");
1985 }
1986
1987 }
1988
1989 void
1990 Data::extractData(const std::string fileName,
1991 const FunctionSpace& fspace)
1992 {
1993 //
1994 // Can only extract Data to an object which is initially DataEmpty
1995 if (!isEmpty()) {
1996 throw DataException("extractData Error: can only extract to DataEmpty object");
1997 }
1998
1999 cout << "Extracting Data object from: " << fileName << endl;
2000
2001 int dataType;
2002 int noSamples;
2003 int noDPPSample;
2004 int functionSpaceType;
2005 int dataPointRank;
2006 int dataPointSize;
2007 int dataLength;
2008 DataArrayView::ShapeType dataPointShape;
2009 int flatShape[4];
2010
2011 //
2012 // Open the archive file
2013 ifstream archiveFile;
2014 archiveFile.open(fileName.data(), ios::in);
2015
2016 if (!archiveFile.good()) {
2017 throw DataException("extractData Error: problem opening archive file");
2018 }
2019
2020 //
2021 // Read common data items from archive file
2022 archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
2023 archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2024 archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2025 archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2026 archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2027 archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2028 archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
2029 for (int dim = 0; dim < 4; dim++) {
2030 archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2031 if (flatShape[dim]>0) {
2032 dataPointShape.push_back(flatShape[dim]);
2033 }
2034 }
2035 int referenceNumbers[noSamples];
2036 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2037 archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2038 }
2039 int tagNumbers[noSamples];
2040 if (dataType==2) {
2041 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2042 archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2043 }
2044 }
2045
2046 if (!archiveFile.good()) {
2047 throw DataException("extractData Error: problem reading from archive file");
2048 }
2049
2050 //
2051 // Verify the values just read from the archive file
2052 switch (dataType) {
2053 case 0:
2054 cout << "\tdataType: DataEmpty" << endl;
2055 break;
2056 case 1:
2057 cout << "\tdataType: DataConstant" << endl;
2058 break;
2059 case 2:
2060 cout << "\tdataType: DataTagged" << endl;
2061 break;
2062 case 3:
2063 cout << "\tdataType: DataExpanded" << endl;
2064 break;
2065 default:
2066 throw DataException("extractData Error: undefined dataType read from archive file");
2067 break;
2068 }
2069
2070 cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2071 cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2072 cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2073 cout << "\tshape: < ";
2074 for (int dim = 0; dim < dataPointRank; dim++) {
2075 cout << dataPointShape[dim] << " ";
2076 }
2077 cout << ">" << endl;
2078
2079 //
2080 // Verify that supplied FunctionSpace object is compatible with this Data object.
2081 if ( (fspace.getTypeCode()!=functionSpaceType) ||
2082 (fspace.getNumSamples()!=noSamples) ||
2083 (fspace.getNumDPPSample()!=noDPPSample)
2084 ) {
2085 throw DataException("extractData Error: incompatible FunctionSpace");
2086 }
2087 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2088 if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2089 throw DataException("extractData Error: incompatible FunctionSpace");
2090 }
2091 }
2092 if (dataType==2) {
2093 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2094 if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2095 throw DataException("extractData Error: incompatible FunctionSpace");
2096 }
2097 }
2098 }
2099
2100 //
2101 // Construct a DataVector to hold underlying data values
2102 DataVector dataVec(dataLength);
2103
2104 //
2105 // Load this DataVector with the appropriate values
2106 int noValues;
2107 archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2108 cout << "\tnoValues: " << noValues << endl;
2109 switch (dataType) {
2110 case 0:
2111 // DataEmpty
2112 if (noValues != 0) {
2113 throw DataException("extractData Error: problem reading data from archive file");
2114 }
2115 break;
2116 case 1:
2117 // DataConstant
2118 if (dataVec.extractData(archiveFile,noValues)) {
2119 throw DataException("extractData Error: problem reading data from archive file");
2120 }
2121 break;
2122 case 2:
2123 // DataTagged
2124 if (dataVec.extractData(archiveFile,noValues)) {
2125 throw DataException("extractData Error: problem reading data from archive file");
2126 }
2127 break;
2128 case 3:
2129 // DataExpanded
2130 if (dataVec.extractData(archiveFile,noValues)) {
2131 throw DataException("extractData Error: problem reading data from archive file");
2132 }
2133 break;
2134 }
2135
2136 if (!archiveFile.good()) {
2137 throw DataException("extractData Error: problem reading from archive file");
2138 }
2139
2140 //
2141 // Close archive file
2142 archiveFile.close();
2143
2144 if (!archiveFile.good()) {
2145 throw DataException("extractData Error: problem closing archive file");
2146 }
2147
2148 //
2149 // Construct an appropriate Data object
2150 DataAbstract* tempData;
2151 switch (dataType) {
2152 case 0:
2153 // DataEmpty
2154 tempData=new DataEmpty();
2155 break;
2156 case 1:
2157 // DataConstant
2158 tempData=new DataConstant(fspace,dataPointShape,dataVec);
2159 break;
2160 case 2:
2161 // DataTagged
2162 tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2163 break;
2164 case 3:
2165 // DataExpanded
2166 tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2167 break;
2168 }
2169 shared_ptr<DataAbstract> temp_data(tempData);
2170 m_data=temp_data;
2171 }
2172
2173 ostream& escript::operator<<(ostream& o, const Data& data)
2174 {
2175 o << data.toString();
2176 return o;
2177 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26