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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26