/[escript]/branches/intelc_win32/escript/src/Data.cpp
ViewVC logotype

Contents of /branches/intelc_win32/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26