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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26