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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26