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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 783 - (show annotations)
Tue Jul 18 01:32:50 2006 UTC (13 years, 6 months ago) by gross
File size: 62557 byte(s)
coordinates, element size and normals returned by corresponding
FunctionSpace mesthods are now protected against updates. So 
+=, -=, *=, /=, setTaggedValue, fillFromNumArray will through an
excpetion.

The FunctionSpace class does nut buffer the oordinates, element size and
normals yet.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26