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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 782 - (show annotations)
Tue Jul 18 00:47:47 2006 UTC (12 years, 9 months ago) by bcumming
File size: 60681 byte(s)
Large number of changes to Finley for meshing in MPI.

- optimisation and neatening up of rectcanglular mesh generation code
- first and second order 1D, 2D and 3D rectangular meshes are now
  available in finley and escript using MPI.
- reduced meshes now generated in MPI, and interpolation to and from 
  reduced data types now supported.  

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26