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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 800 - (show annotations)
Tue Aug 8 11:23:18 2006 UTC (12 years, 10 months ago) by gross
File size: 65317 byte(s)
new function _swap. Python wrapper + testing is still missing.


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 Data::convertToNumArrayFromDPNo(int procNo,
786 int sampleNo,
787 int dataPointNo)
788
789 {
790 size_t length=0;
791 int i, j, k, l, pos;
792
793 //
794 // Check a valid sample number has been supplied
795 if (sampleNo >= getNumSamples()) {
796 throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
797 }
798
799 //
800 // Check a valid data point number has been supplied
801 if (dataPointNo >= getNumDataPointsPerSample()) {
802 throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");
803 }
804
805 //
806 // determine the rank and shape of each data point
807 int dataPointRank = getDataPointRank();
808 DataArrayView::ShapeType dataPointShape = getDataPointShape();
809
810 //
811 // create the numeric array to be returned
812 boost::python::numeric::array numArray(0.0);
813
814 //
815 // the shape of the returned numeric array will be the same
816 // as that of the data point
817 int arrayRank = dataPointRank;
818 DataArrayView::ShapeType arrayShape = dataPointShape;
819
820 //
821 // resize the numeric array to the shape just calculated
822 if (arrayRank==0) {
823 numArray.resize(1);
824 }
825 if (arrayRank==1) {
826 numArray.resize(arrayShape[0]);
827 }
828 if (arrayRank==2) {
829 numArray.resize(arrayShape[0],arrayShape[1]);
830 }
831 if (arrayRank==3) {
832 numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
833 }
834 if (arrayRank==4) {
835 numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
836 }
837
838 // added for the MPI communication
839 length=1;
840 for( i=0; i<arrayRank; i++ )
841 length *= arrayShape[i];
842 double *tmpData = new double[length];
843
844 //
845 // load the values for the data point into the numeric array.
846
847 // updated for the MPI case
848 if( get_MPIRank()==procNo ){
849 // create a view of the data if it is stored locally
850 DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
851
852 // pack the data from the view into tmpData for MPI communication
853 pos=0;
854 switch( dataPointRank ){
855 case 0 :
856 tmpData[0] = dataPointView();
857 break;
858 case 1 :
859 for( i=0; i<dataPointShape[0]; i++ )
860 tmpData[i]=dataPointView(i);
861 break;
862 case 2 :
863 for( i=0; i<dataPointShape[0]; i++ )
864 for( j=0; j<dataPointShape[1]; j++, pos++ )
865 tmpData[pos]=dataPointView(i,j);
866 break;
867 case 3 :
868 for( i=0; i<dataPointShape[0]; i++ )
869 for( j=0; j<dataPointShape[1]; j++ )
870 for( k=0; k<dataPointShape[2]; k++, pos++ )
871 tmpData[pos]=dataPointView(i,j,k);
872 break;
873 case 4 :
874 for( i=0; i<dataPointShape[0]; i++ )
875 for( j=0; j<dataPointShape[1]; j++ )
876 for( k=0; k<dataPointShape[2]; k++ )
877 for( l=0; l<dataPointShape[3]; l++, pos++ )
878 tmpData[pos]=dataPointView(i,j,k,l);
879 break;
880 }
881 }
882 #ifdef PASO_MPI
883 // broadcast the data to all other processes
884 MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
885 #endif
886
887 // unpack the data
888 switch( dataPointRank ){
889 case 0 :
890 numArray[i]=tmpData[0];
891 break;
892 case 1 :
893 for( i=0; i<dataPointShape[0]; i++ )
894 numArray[i]=tmpData[i];
895 break;
896 case 2 :
897 for( i=0; i<dataPointShape[0]; i++ )
898 for( j=0; j<dataPointShape[1]; j++ )
899 tmpData[i+j*dataPointShape[0]];
900 break;
901 case 3 :
902 for( i=0; i<dataPointShape[0]; i++ )
903 for( j=0; j<dataPointShape[1]; j++ )
904 for( k=0; k<dataPointShape[2]; k++ )
905 tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
906 break;
907 case 4 :
908 for( i=0; i<dataPointShape[0]; i++ )
909 for( j=0; j<dataPointShape[1]; j++ )
910 for( k=0; k<dataPointShape[2]; k++ )
911 for( l=0; l<dataPointShape[3]; l++ )
912 tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
913 break;
914 }
915
916 delete [] tmpData;
917 /*
918 if (dataPointRank==0) {
919 numArray[0]=dataPointView();
920 }
921 if (dataPointRank==1) {
922 for (int i=0; i<dataPointShape[0]; i++) {
923 numArray[i]=dataPointView(i);
924 }
925 }
926 if (dataPointRank==2) {
927 for (int i=0; i<dataPointShape[0]; i++) {
928 for (int j=0; j<dataPointShape[1]; j++) {
929 numArray[i][j] = dataPointView(i,j);
930 }
931 }
932 }
933 if (dataPointRank==3) {
934 for (int i=0; i<dataPointShape[0]; i++) {
935 for (int j=0; j<dataPointShape[1]; j++) {
936 for (int k=0; k<dataPointShape[2]; k++) {
937 numArray[i][j][k]=dataPointView(i,j,k);
938 }
939 }
940 }
941 }
942 if (dataPointRank==4) {
943 for (int i=0; i<dataPointShape[0]; i++) {
944 for (int j=0; j<dataPointShape[1]; j++) {
945 for (int k=0; k<dataPointShape[2]; k++) {
946 for (int l=0; l<dataPointShape[3]; l++) {
947 numArray[i][j][k][l]=dataPointView(i,j,k,l);
948 }
949 }
950 }
951 }
952 }
953 */
954
955 //
956 // return the loaded array
957 return numArray;
958 }
959
960 boost::python::numeric::array
961 Data::integrate() const
962 {
963 int index;
964 int rank = getDataPointRank();
965 DataArrayView::ShapeType shape = getDataPointShape();
966
967 #if defined DOPROF
968 profData->integrate++;
969 #endif
970
971 //
972 // calculate the integral values
973 vector<double> integrals(getDataPointSize());
974 AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);
975
976 //
977 // create the numeric array to be returned
978 // and load the array with the integral values
979 boost::python::numeric::array bp_array(1.0);
980 if (rank==0) {
981 bp_array.resize(1);
982 index = 0;
983 bp_array[0] = integrals[index];
984 }
985 if (rank==1) {
986 bp_array.resize(shape[0]);
987 for (int i=0; i<shape[0]; i++) {
988 index = i;
989 bp_array[i] = integrals[index];
990 }
991 }
992 if (rank==2) {
993 bp_array.resize(shape[0],shape[1]);
994 for (int i=0; i<shape[0]; i++) {
995 for (int j=0; j<shape[1]; j++) {
996 index = i + shape[0] * j;
997 bp_array[make_tuple(i,j)] = integrals[index];
998 }
999 }
1000 }
1001 if (rank==3) {
1002 bp_array.resize(shape[0],shape[1],shape[2]);
1003 for (int i=0; i<shape[0]; i++) {
1004 for (int j=0; j<shape[1]; j++) {
1005 for (int k=0; k<shape[2]; k++) {
1006 index = i + shape[0] * ( j + shape[1] * k );
1007 bp_array[make_tuple(i,j,k)] = integrals[index];
1008 }
1009 }
1010 }
1011 }
1012 if (rank==4) {
1013 bp_array.resize(shape[0],shape[1],shape[2],shape[3]);
1014 for (int i=0; i<shape[0]; i++) {
1015 for (int j=0; j<shape[1]; j++) {
1016 for (int k=0; k<shape[2]; k++) {
1017 for (int l=0; l<shape[3]; l++) {
1018 index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
1019 bp_array[make_tuple(i,j,k,l)] = integrals[index];
1020 }
1021 }
1022 }
1023 }
1024 }
1025
1026 //
1027 // return the loaded array
1028 return bp_array;
1029 }
1030
1031 Data
1032 Data::sin() const
1033 {
1034 #if defined DOPROF
1035 profData->unary++;
1036 #endif
1037 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
1038 }
1039
1040 Data
1041 Data::cos() const
1042 {
1043 #if defined DOPROF
1044 profData->unary++;
1045 #endif
1046 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
1047 }
1048
1049 Data
1050 Data::tan() const
1051 {
1052 #if defined DOPROF
1053 profData->unary++;
1054 #endif
1055 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
1056 }
1057
1058 Data
1059 Data::asin() const
1060 {
1061 #if defined DOPROF
1062 profData->unary++;
1063 #endif
1064 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
1065 }
1066
1067 Data
1068 Data::acos() const
1069 {
1070 #if defined DOPROF
1071 profData->unary++;
1072 #endif
1073 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
1074 }
1075
1076 Data
1077 Data::atan() const
1078 {
1079 #if defined DOPROF
1080 profData->unary++;
1081 #endif
1082 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
1083 }
1084
1085 Data
1086 Data::sinh() const
1087 {
1088 #if defined DOPROF
1089 profData->unary++;
1090 #endif
1091 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
1092 }
1093
1094 Data
1095 Data::cosh() const
1096 {
1097 #if defined DOPROF
1098 profData->unary++;
1099 #endif
1100 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
1101 }
1102
1103 Data
1104 Data::tanh() const
1105 {
1106 #if defined DOPROF
1107 profData->unary++;
1108 #endif
1109 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
1110 }
1111
1112 Data
1113 Data::asinh() const
1114 {
1115 #if defined DOPROF
1116 profData->unary++;
1117 #endif
1118 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
1119 }
1120
1121 Data
1122 Data::acosh() const
1123 {
1124 #if defined DOPROF
1125 profData->unary++;
1126 #endif
1127 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
1128 }
1129
1130 Data
1131 Data::atanh() const
1132 {
1133 #if defined DOPROF
1134 profData->unary++;
1135 #endif
1136 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
1137 }
1138
1139 Data
1140 Data::log10() const
1141 {
1142 #if defined DOPROF
1143 profData->unary++;
1144 #endif
1145 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
1146 }
1147
1148 Data
1149 Data::log() const
1150 {
1151 #if defined DOPROF
1152 profData->unary++;
1153 #endif
1154 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
1155 }
1156
1157 Data
1158 Data::sign() const
1159 {
1160 #if defined DOPROF
1161 profData->unary++;
1162 #endif
1163 return escript::unaryOp(*this,escript::fsign);
1164 }
1165
1166 Data
1167 Data::abs() const
1168 {
1169 #if defined DOPROF
1170 profData->unary++;
1171 #endif
1172 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
1173 }
1174
1175 Data
1176 Data::neg() const
1177 {
1178 #if defined DOPROF
1179 profData->unary++;
1180 #endif
1181 return escript::unaryOp(*this,negate<double>());
1182 }
1183
1184 Data
1185 Data::pos() const
1186 {
1187 #if defined DOPROF
1188 profData->unary++;
1189 #endif
1190 Data result;
1191 // perform a deep copy
1192 result.copy(*this);
1193 return result;
1194 }
1195
1196 Data
1197 Data::exp() const
1198 {
1199 #if defined DOPROF
1200 profData->unary++;
1201 #endif
1202 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1203 }
1204
1205 Data
1206 Data::sqrt() const
1207 {
1208 #if defined DOPROF
1209 profData->unary++;
1210 #endif
1211 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1212 }
1213
1214 double
1215 Data::Lsup() const
1216 {
1217 double localValue, globalValue;
1218 #if defined DOPROF
1219 profData->reduction1++;
1220 #endif
1221 //
1222 // set the initial absolute maximum value to zero
1223
1224 AbsMax abs_max_func;
1225 localValue = algorithm(abs_max_func,0);
1226 #ifdef PASO_MPI
1227 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1228 return globalValue;
1229 #else
1230 return localValue;
1231 #endif
1232 }
1233
1234 double
1235 Data::Linf() const
1236 {
1237 double localValue, globalValue;
1238 #if defined DOPROF
1239 profData->reduction1++;
1240 #endif
1241 //
1242 // set the initial absolute minimum value to max double
1243 AbsMin abs_min_func;
1244 localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1245
1246 #ifdef PASO_MPI
1247 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1248 return globalValue;
1249 #else
1250 return localValue;
1251 #endif
1252 }
1253
1254 double
1255 Data::sup() const
1256 {
1257 double localValue, globalValue;
1258 #if defined DOPROF
1259 profData->reduction1++;
1260 #endif
1261 //
1262 // set the initial maximum value to min possible double
1263 FMax fmax_func;
1264 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1265 #ifdef PASO_MPI
1266 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1267 return globalValue;
1268 #else
1269 return localValue;
1270 #endif
1271 }
1272
1273 double
1274 Data::inf() const
1275 {
1276 double localValue, globalValue;
1277 #if defined DOPROF
1278 profData->reduction1++;
1279 #endif
1280 //
1281 // set the initial minimum value to max possible double
1282 FMin fmin_func;
1283 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1284 #ifdef PASO_MPI
1285 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1286 return globalValue;
1287 #else
1288 return localValue;
1289 #endif
1290 }
1291
1292 /* TODO */
1293 /* global reduction */
1294 Data
1295 Data::maxval() const
1296 {
1297 #if defined DOPROF
1298 profData->reduction2++;
1299 #endif
1300 //
1301 // set the initial maximum value to min possible double
1302 FMax fmax_func;
1303 return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1304 }
1305
1306 Data
1307 Data::minval() const
1308 {
1309 #if defined DOPROF
1310 profData->reduction2++;
1311 #endif
1312 //
1313 // set the initial minimum value to max possible double
1314 FMin fmin_func;
1315 return dp_algorithm(fmin_func,numeric_limits<double>::max());
1316 }
1317
1318 Data
1319 Data::swap(int axis_offset) const
1320 {
1321 #if defined DOPROF
1322 profData->unary++;
1323 #endif
1324 DataArrayView::ShapeType s=getDataPointShape();
1325 DataArrayView::ShapeType ev_shape;
1326 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1327 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1328 int rank=getDataPointRank();
1329 if (axis_offset<0 || axis_offset+1>rank) {
1330 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank-1=" + rank-1);
1331 }
1332 for (int i=0; i<rank; i++) {
1333 if (i == axis_offset) {
1334 ev_shape.push_back(s[axis_offset+1]);
1335 } else if (i == axis_offset+1) {
1336 ev_shape.push_back(s[axis_offset]);
1337 } else {
1338 ev_shape.push_back(s[i]);
1339 }
1340 }
1341 Data ev(0.,ev_shape,getFunctionSpace());
1342 ev.typeMatchRight(*this);
1343 m_data->swap(ev.m_data.get(), axis_offset);
1344 return ev;
1345
1346 }
1347
1348 Data
1349 Data::symmetric() const
1350 {
1351 #if defined DOPROF
1352 profData->unary++;
1353 #endif
1354 // check input
1355 DataArrayView::ShapeType s=getDataPointShape();
1356 if (getDataPointRank()==2) {
1357 if(s[0] != s[1])
1358 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1359 }
1360 else if (getDataPointRank()==4) {
1361 if(!(s[0] == s[2] && s[1] == s[3]))
1362 throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1363 }
1364 else {
1365 throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1366 }
1367 Data ev(0.,getDataPointShape(),getFunctionSpace());
1368 ev.typeMatchRight(*this);
1369 m_data->symmetric(ev.m_data.get());
1370 return ev;
1371 }
1372
1373 Data
1374 Data::nonsymmetric() const
1375 {
1376 #if defined DOPROF
1377 profData->unary++;
1378 #endif
1379 // check input
1380 DataArrayView::ShapeType s=getDataPointShape();
1381 if (getDataPointRank()==2) {
1382 if(s[0] != s[1])
1383 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1384 DataArrayView::ShapeType ev_shape;
1385 ev_shape.push_back(s[0]);
1386 ev_shape.push_back(s[1]);
1387 Data ev(0.,ev_shape,getFunctionSpace());
1388 ev.typeMatchRight(*this);
1389 m_data->nonsymmetric(ev.m_data.get());
1390 return ev;
1391 }
1392 else if (getDataPointRank()==4) {
1393 if(!(s[0] == s[2] && s[1] == s[3]))
1394 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1395 DataArrayView::ShapeType ev_shape;
1396 ev_shape.push_back(s[0]);
1397 ev_shape.push_back(s[1]);
1398 ev_shape.push_back(s[2]);
1399 ev_shape.push_back(s[3]);
1400 Data ev(0.,ev_shape,getFunctionSpace());
1401 ev.typeMatchRight(*this);
1402 m_data->nonsymmetric(ev.m_data.get());
1403 return ev;
1404 }
1405 else {
1406 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1407 }
1408 }
1409
1410 Data
1411 Data::trace(int axis_offset) const
1412 {
1413 #if defined DOPROF
1414 profData->unary++;
1415 #endif
1416 DataArrayView::ShapeType s=getDataPointShape();
1417 if (getDataPointRank()==2) {
1418 DataArrayView::ShapeType ev_shape;
1419 Data ev(0.,ev_shape,getFunctionSpace());
1420 ev.typeMatchRight(*this);
1421 m_data->trace(ev.m_data.get(), axis_offset);
1422 return ev;
1423 }
1424 if (getDataPointRank()==3) {
1425 DataArrayView::ShapeType ev_shape;
1426 if (axis_offset==0) {
1427 int s2=s[2];
1428 ev_shape.push_back(s2);
1429 }
1430 else if (axis_offset==1) {
1431 int s0=s[0];
1432 ev_shape.push_back(s0);
1433 }
1434 Data ev(0.,ev_shape,getFunctionSpace());
1435 ev.typeMatchRight(*this);
1436 m_data->trace(ev.m_data.get(), axis_offset);
1437 return ev;
1438 }
1439 if (getDataPointRank()==4) {
1440 DataArrayView::ShapeType ev_shape;
1441 if (axis_offset==0) {
1442 ev_shape.push_back(s[2]);
1443 ev_shape.push_back(s[3]);
1444 }
1445 else if (axis_offset==1) {
1446 ev_shape.push_back(s[0]);
1447 ev_shape.push_back(s[3]);
1448 }
1449 else if (axis_offset==2) {
1450 ev_shape.push_back(s[0]);
1451 ev_shape.push_back(s[1]);
1452 }
1453 Data ev(0.,ev_shape,getFunctionSpace());
1454 ev.typeMatchRight(*this);
1455 m_data->trace(ev.m_data.get(), axis_offset);
1456 return ev;
1457 }
1458 else {
1459 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1460 }
1461 }
1462
1463 Data
1464 Data::transpose(int axis_offset) const
1465 {
1466 #if defined DOPROF
1467 profData->unary++;
1468 #endif
1469 DataArrayView::ShapeType s=getDataPointShape();
1470 DataArrayView::ShapeType ev_shape;
1471 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1472 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1473 int rank=getDataPointRank();
1474 if (axis_offset<0 || axis_offset>rank) {
1475 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1476 }
1477 for (int i=0; i<rank; i++) {
1478 int index = (axis_offset+i)%rank;
1479 ev_shape.push_back(s[index]); // Append to new shape
1480 }
1481 Data ev(0.,ev_shape,getFunctionSpace());
1482 ev.typeMatchRight(*this);
1483 m_data->transpose(ev.m_data.get(), axis_offset);
1484 return ev;
1485 }
1486
1487 Data
1488 Data::eigenvalues() const
1489 {
1490 #if defined DOPROF
1491 profData->unary++;
1492 #endif
1493 // check input
1494 DataArrayView::ShapeType s=getDataPointShape();
1495 if (getDataPointRank()!=2)
1496 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1497 if(s[0] != s[1])
1498 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1499 // create return
1500 DataArrayView::ShapeType ev_shape(1,s[0]);
1501 Data ev(0.,ev_shape,getFunctionSpace());
1502 ev.typeMatchRight(*this);
1503 m_data->eigenvalues(ev.m_data.get());
1504 return ev;
1505 }
1506
1507 const boost::python::tuple
1508 Data::eigenvalues_and_eigenvectors(const double tol) const
1509 {
1510 #if defined DOPROF
1511 profData->unary++;
1512 #endif
1513 DataArrayView::ShapeType s=getDataPointShape();
1514 if (getDataPointRank()!=2)
1515 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1516 if(s[0] != s[1])
1517 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1518 // create return
1519 DataArrayView::ShapeType ev_shape(1,s[0]);
1520 Data ev(0.,ev_shape,getFunctionSpace());
1521 ev.typeMatchRight(*this);
1522 DataArrayView::ShapeType V_shape(2,s[0]);
1523 Data V(0.,V_shape,getFunctionSpace());
1524 V.typeMatchRight(*this);
1525 m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1526 return make_tuple(boost::python::object(ev),boost::python::object(V));
1527 }
1528
1529 const boost::python::tuple
1530 Data::mindp() const
1531 {
1532 // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1533 // abort (for unknown reasons) if there are openmp directives with it in the
1534 // surrounding function
1535
1536 int SampleNo;
1537 int DataPointNo;
1538 int ProcNo;
1539 calc_mindp(ProcNo,SampleNo,DataPointNo);
1540 return make_tuple(ProcNo,SampleNo,DataPointNo);
1541 }
1542
1543 void
1544 Data::calc_mindp( int& ProcNo,
1545 int& SampleNo,
1546 int& DataPointNo) const
1547 {
1548 int i,j;
1549 int lowi=0,lowj=0;
1550 double min=numeric_limits<double>::max();
1551
1552 Data temp=minval();
1553
1554 int numSamples=temp.getNumSamples();
1555 int numDPPSample=temp.getNumDataPointsPerSample();
1556
1557 double next,local_min;
1558 int local_lowi,local_lowj;
1559
1560 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1561 {
1562 local_min=min;
1563 #pragma omp for private(i,j) schedule(static)
1564 for (i=0; i<numSamples; i++) {
1565 for (j=0; j<numDPPSample; j++) {
1566 next=temp.getDataPoint(i,j)();
1567 if (next<local_min) {
1568 local_min=next;
1569 local_lowi=i;
1570 local_lowj=j;
1571 }
1572 }
1573 }
1574 #pragma omp critical
1575 if (local_min<min) {
1576 min=local_min;
1577 lowi=local_lowi;
1578 lowj=local_lowj;
1579 }
1580 }
1581
1582 #ifdef PASO_MPI
1583 // determine the processor on which the minimum occurs
1584 next = temp.getDataPoint(lowi,lowj)();
1585 int lowProc = 0;
1586 double *globalMins = new double[get_MPISize()+1];
1587 int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1588
1589 if( get_MPIRank()==0 ){
1590 next = globalMins[lowProc];
1591 for( i=1; i<get_MPISize(); i++ )
1592 if( next>globalMins[i] ){
1593 lowProc = i;
1594 next = globalMins[i];
1595 }
1596 }
1597 MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1598
1599 delete [] globalMins;
1600 ProcNo = lowProc;
1601 #else
1602 ProcNo = 0;
1603 #endif
1604 SampleNo = lowi;
1605 DataPointNo = lowj;
1606 }
1607
1608 void
1609 Data::saveDX(std::string fileName) const
1610 {
1611 boost::python::dict args;
1612 args["data"]=boost::python::object(this);
1613 getDomain().saveDX(fileName,args);
1614 return;
1615 }
1616
1617 void
1618 Data::saveVTK(std::string fileName) const
1619 {
1620 boost::python::dict args;
1621 args["data"]=boost::python::object(this);
1622 getDomain().saveVTK(fileName,args);
1623 return;
1624 }
1625
1626 Data&
1627 Data::operator+=(const Data& right)
1628 {
1629 if (isProtected()) {
1630 throw DataException("Error - attempt to update protected Data object.");
1631 }
1632 #if defined DOPROF
1633 profData->binary++;
1634 #endif
1635 binaryOp(right,plus<double>());
1636 return (*this);
1637 }
1638
1639 Data&
1640 Data::operator+=(const boost::python::object& right)
1641 {
1642 if (isProtected()) {
1643 throw DataException("Error - attempt to update protected Data object.");
1644 }
1645 #if defined DOPROF
1646 profData->binary++;
1647 #endif
1648 binaryOp(right,plus<double>());
1649 return (*this);
1650 }
1651
1652 Data&
1653 Data::operator-=(const Data& right)
1654 {
1655 if (isProtected()) {
1656 throw DataException("Error - attempt to update protected Data object.");
1657 }
1658 #if defined DOPROF
1659 profData->binary++;
1660 #endif
1661 binaryOp(right,minus<double>());
1662 return (*this);
1663 }
1664
1665 Data&
1666 Data::operator-=(const boost::python::object& right)
1667 {
1668 if (isProtected()) {
1669 throw DataException("Error - attempt to update protected Data object.");
1670 }
1671 #if defined DOPROF
1672 profData->binary++;
1673 #endif
1674 binaryOp(right,minus<double>());
1675 return (*this);
1676 }
1677
1678 Data&
1679 Data::operator*=(const Data& right)
1680 {
1681 if (isProtected()) {
1682 throw DataException("Error - attempt to update protected Data object.");
1683 }
1684 #if defined DOPROF
1685 profData->binary++;
1686 #endif
1687 binaryOp(right,multiplies<double>());
1688 return (*this);
1689 }
1690
1691 Data&
1692 Data::operator*=(const boost::python::object& right)
1693 {
1694 if (isProtected()) {
1695 throw DataException("Error - attempt to update protected Data object.");
1696 }
1697 #if defined DOPROF
1698 profData->binary++;
1699 #endif
1700 binaryOp(right,multiplies<double>());
1701 return (*this);
1702 }
1703
1704 Data&
1705 Data::operator/=(const Data& right)
1706 {
1707 if (isProtected()) {
1708 throw DataException("Error - attempt to update protected Data object.");
1709 }
1710 #if defined DOPROF
1711 profData->binary++;
1712 #endif
1713 binaryOp(right,divides<double>());
1714 return (*this);
1715 }
1716
1717 Data&
1718 Data::operator/=(const boost::python::object& right)
1719 {
1720 if (isProtected()) {
1721 throw DataException("Error - attempt to update protected Data object.");
1722 }
1723 #if defined DOPROF
1724 profData->binary++;
1725 #endif
1726 binaryOp(right,divides<double>());
1727 return (*this);
1728 }
1729
1730 Data
1731 Data::rpowO(const boost::python::object& left) const
1732 {
1733 if (isProtected()) {
1734 throw DataException("Error - attempt to update protected Data object.");
1735 }
1736 #if defined DOPROF
1737 profData->binary++;
1738 #endif
1739 Data left_d(left,*this);
1740 return left_d.powD(*this);
1741 }
1742
1743 Data
1744 Data::powO(const boost::python::object& right) const
1745 {
1746 #if defined DOPROF
1747 profData->binary++;
1748 #endif
1749 Data result;
1750 result.copy(*this);
1751 result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1752 return result;
1753 }
1754
1755 Data
1756 Data::powD(const Data& right) const
1757 {
1758 #if defined DOPROF
1759 profData->binary++;
1760 #endif
1761 Data result;
1762 result.copy(*this);
1763 result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1764 return result;
1765 }
1766
1767
1768 //
1769 // NOTE: It is essential to specify the namespace this operator belongs to
1770 Data
1771 escript::operator+(const Data& left, const Data& right)
1772 {
1773 Data result;
1774 //
1775 // perform a deep copy
1776 result.copy(left);
1777 result+=right;
1778 return result;
1779 }
1780
1781 //
1782 // NOTE: It is essential to specify the namespace this operator belongs to
1783 Data
1784 escript::operator-(const Data& left, const Data& right)
1785 {
1786 Data result;
1787 //
1788 // perform a deep copy
1789 result.copy(left);
1790 result-=right;
1791 return result;
1792 }
1793
1794 //
1795 // NOTE: It is essential to specify the namespace this operator belongs to
1796 Data
1797 escript::operator*(const Data& left, const Data& right)
1798 {
1799 Data result;
1800 //
1801 // perform a deep copy
1802 result.copy(left);
1803 result*=right;
1804 return result;
1805 }
1806
1807 //
1808 // NOTE: It is essential to specify the namespace this operator belongs to
1809 Data
1810 escript::operator/(const Data& left, const Data& right)
1811 {
1812 Data result;
1813 //
1814 // perform a deep copy
1815 result.copy(left);
1816 result/=right;
1817 return result;
1818 }
1819
1820 //
1821 // NOTE: It is essential to specify the namespace this operator belongs to
1822 Data
1823 escript::operator+(const Data& left, const boost::python::object& right)
1824 {
1825 //
1826 // Convert to DataArray format if possible
1827 DataArray temp(right);
1828 Data result;
1829 //
1830 // perform a deep copy
1831 result.copy(left);
1832 result+=right;
1833 return result;
1834 }
1835
1836 //
1837 // NOTE: It is essential to specify the namespace this operator belongs to
1838 Data
1839 escript::operator-(const Data& left, const boost::python::object& right)
1840 {
1841 //
1842 // Convert to DataArray format if possible
1843 DataArray temp(right);
1844 Data result;
1845 //
1846 // perform a deep copy
1847 result.copy(left);
1848 result-=right;
1849 return result;
1850 }
1851
1852 //
1853 // NOTE: It is essential to specify the namespace this operator belongs to
1854 Data
1855 escript::operator*(const Data& left, const boost::python::object& right)
1856 {
1857 //
1858 // Convert to DataArray format if possible
1859 DataArray temp(right);
1860 Data result;
1861 //
1862 // perform a deep copy
1863 result.copy(left);
1864 result*=right;
1865 return result;
1866 }
1867
1868 //
1869 // NOTE: It is essential to specify the namespace this operator belongs to
1870 Data
1871 escript::operator/(const Data& left, const boost::python::object& right)
1872 {
1873 //
1874 // Convert to DataArray format if possible
1875 DataArray temp(right);
1876 Data result;
1877 //
1878 // perform a deep copy
1879 result.copy(left);
1880 result/=right;
1881 return result;
1882 }
1883
1884 //
1885 // NOTE: It is essential to specify the namespace this operator belongs to
1886 Data
1887 escript::operator+(const boost::python::object& left, const Data& right)
1888 {
1889 //
1890 // Construct the result using the given value and the other parameters
1891 // from right
1892 Data result(left,right);
1893 result+=right;
1894 return result;
1895 }
1896
1897 //
1898 // NOTE: It is essential to specify the namespace this operator belongs to
1899 Data
1900 escript::operator-(const boost::python::object& left, const Data& right)
1901 {
1902 //
1903 // Construct the result using the given value and the other parameters
1904 // from right
1905 Data result(left,right);
1906 result-=right;
1907 return result;
1908 }
1909
1910 //
1911 // NOTE: It is essential to specify the namespace this operator belongs to
1912 Data
1913 escript::operator*(const boost::python::object& left, const Data& right)
1914 {
1915 //
1916 // Construct the result using the given value and the other parameters
1917 // from right
1918 Data result(left,right);
1919 result*=right;
1920 return result;
1921 }
1922
1923 //
1924 // NOTE: It is essential to specify the namespace this operator belongs to
1925 Data
1926 escript::operator/(const boost::python::object& left, const Data& right)
1927 {
1928 //
1929 // Construct the result using the given value and the other parameters
1930 // from right
1931 Data result(left,right);
1932 result/=right;
1933 return result;
1934 }
1935
1936 //
1937 //bool escript::operator==(const Data& left, const Data& right)
1938 //{
1939 // /*
1940 // NB: this operator does very little at this point, and isn't to
1941 // be relied on. Requires further implementation.
1942 // */
1943 //
1944 // bool ret;
1945 //
1946 // if (left.isEmpty()) {
1947 // if(!right.isEmpty()) {
1948 // ret = false;
1949 // } else {
1950 // ret = true;
1951 // }
1952 // }
1953 //
1954 // if (left.isConstant()) {
1955 // if(!right.isConstant()) {
1956 // ret = false;
1957 // } else {
1958 // ret = true;
1959 // }
1960 // }
1961 //
1962 // if (left.isTagged()) {
1963 // if(!right.isTagged()) {
1964 // ret = false;
1965 // } else {
1966 // ret = true;
1967 // }
1968 // }
1969 //
1970 // if (left.isExpanded()) {
1971 // if(!right.isExpanded()) {
1972 // ret = false;
1973 // } else {
1974 // ret = true;
1975 // }
1976 // }
1977 //
1978 // return ret;
1979 //}
1980
1981 /* TODO */
1982 /* global reduction */
1983 Data
1984 Data::getItem(const boost::python::object& key) const
1985 {
1986 const DataArrayView& view=getPointDataView();
1987
1988 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1989
1990 if (slice_region.size()!=view.getRank()) {
1991 throw DataException("Error - slice size does not match Data rank.");
1992 }
1993
1994 return getSlice(slice_region);
1995 }
1996
1997 /* TODO */
1998 /* global reduction */
1999 Data
2000 Data::getSlice(const DataArrayView::RegionType& region) const
2001 {
2002 #if defined DOPROF
2003 profData->slicing++;
2004 #endif
2005 return Data(*this,region);
2006 }
2007
2008 /* TODO */
2009 /* global reduction */
2010 void
2011 Data::setItemO(const boost::python::object& key,
2012 const boost::python::object& value)
2013 {
2014 Data tempData(value,getFunctionSpace());
2015 setItemD(key,tempData);
2016 }
2017
2018 /* TODO */
2019 /* global reduction */
2020 void
2021 Data::setItemD(const boost::python::object& key,
2022 const Data& value)
2023 {
2024 const DataArrayView& view=getPointDataView();
2025
2026 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
2027 if (slice_region.size()!=view.getRank()) {
2028 throw DataException("Error - slice size does not match Data rank.");
2029 }
2030 if (getFunctionSpace()!=value.getFunctionSpace()) {
2031 setSlice(Data(value,getFunctionSpace()),slice_region);
2032 } else {
2033 setSlice(value,slice_region);
2034 }
2035 }
2036
2037 /* TODO */
2038 /* global reduction */
2039 void
2040 Data::setSlice(const Data& value,
2041 const DataArrayView::RegionType& region)
2042 {
2043 if (isProtected()) {
2044 throw DataException("Error - attempt to update protected Data object.");
2045 }
2046 #if defined DOPROF
2047 profData->slicing++;
2048 #endif
2049 Data tempValue(value);
2050 typeMatchLeft(tempValue);
2051 typeMatchRight(tempValue);
2052 m_data->setSlice(tempValue.m_data.get(),region);
2053 }
2054
2055 void
2056 Data::typeMatchLeft(Data& right) const
2057 {
2058 if (isExpanded()){
2059 right.expand();
2060 } else if (isTagged()) {
2061 if (right.isConstant()) {
2062 right.tag();
2063 }
2064 }
2065 }
2066
2067 void
2068 Data::typeMatchRight(const Data& right)
2069 {
2070 if (isTagged()) {
2071 if (right.isExpanded()) {
2072 expand();
2073 }
2074 } else if (isConstant()) {
2075 if (right.isExpanded()) {
2076 expand();
2077 } else if (right.isTagged()) {
2078 tag();
2079 }
2080 }
2081 }
2082
2083 /* TODO */
2084 /* global reduction */
2085 void
2086 Data::setTaggedValue(int tagKey,
2087 const boost::python::object& value)
2088 {
2089 if (isProtected()) {
2090 throw DataException("Error - attempt to update protected Data object.");
2091 }
2092 //
2093 // Ensure underlying data object is of type DataTagged
2094 tag();
2095
2096 if (!isTagged()) {
2097 throw DataException("Error - DataTagged conversion failed!!");
2098 }
2099
2100 //
2101 // Construct DataArray from boost::python::object input value
2102 DataArray valueDataArray(value);
2103
2104 //
2105 // Call DataAbstract::setTaggedValue
2106 m_data->setTaggedValue(tagKey,valueDataArray.getView());
2107 }
2108
2109 /* TODO */
2110 /* global reduction */
2111 void
2112 Data::setTaggedValueFromCPP(int tagKey,
2113 const DataArrayView& value)
2114 {
2115 if (isProtected()) {
2116 throw DataException("Error - attempt to update protected Data object.");
2117 }
2118 //
2119 // Ensure underlying data object is of type DataTagged
2120 tag();
2121
2122 if (!isTagged()) {
2123 throw DataException("Error - DataTagged conversion failed!!");
2124 }
2125
2126 //
2127 // Call DataAbstract::setTaggedValue
2128 m_data->setTaggedValue(tagKey,value);
2129 }
2130
2131 /* TODO */
2132 /* global reduction */
2133 int
2134 Data::getTagNumber(int dpno)
2135 {
2136 return m_data->getTagNumber(dpno);
2137 }
2138
2139 /* TODO */
2140 /* global reduction */
2141 void
2142 Data::setRefValue(int ref,
2143 const boost::python::numeric::array& value)
2144 {
2145 if (isProtected()) {
2146 throw DataException("Error - attempt to update protected Data object.");
2147 }
2148 //
2149 // Construct DataArray from boost::python::object input value
2150 DataArray valueDataArray(value);
2151
2152 //
2153 // Call DataAbstract::setRefValue
2154 m_data->setRefValue(ref,valueDataArray);
2155 }
2156
2157 /* TODO */
2158 /* global reduction */
2159 void
2160 Data::getRefValue(int ref,
2161 boost::python::numeric::array& value)
2162 {
2163 //
2164 // Construct DataArray for boost::python::object return value
2165 DataArray valueDataArray(value);
2166
2167 //
2168 // Load DataArray with values from data-points specified by ref
2169 m_data->getRefValue(ref,valueDataArray);
2170
2171 //
2172 // Load values from valueDataArray into return numarray
2173
2174 // extract the shape of the numarray
2175 int rank = value.getrank();
2176 DataArrayView::ShapeType shape;
2177 for (int i=0; i < rank; i++) {
2178 shape.push_back(extract<int>(value.getshape()[i]));
2179 }
2180
2181 // and load the numarray with the data from the DataArray
2182 DataArrayView valueView = valueDataArray.getView();
2183
2184 if (rank==0) {
2185 boost::python::numeric::array temp_numArray(valueView());
2186 value = temp_numArray;
2187 }
2188 if (rank==1) {
2189 for (int i=0; i < shape[0]; i++) {
2190 value[i] = valueView(i);
2191 }
2192 }
2193 if (rank==2) {
2194 for (int i=0; i < shape[0]; i++) {
2195 for (int j=0; j < shape[1]; j++) {
2196 value[i][j] = valueView(i,j);
2197 }
2198 }
2199 }
2200 if (rank==3) {
2201 for (int i=0; i < shape[0]; i++) {
2202 for (int j=0; j < shape[1]; j++) {
2203 for (int k=0; k < shape[2]; k++) {
2204 value[i][j][k] = valueView(i,j,k);
2205 }
2206 }
2207 }
2208 }
2209 if (rank==4) {
2210 for (int i=0; i < shape[0]; i++) {
2211 for (int j=0; j < shape[1]; j++) {
2212 for (int k=0; k < shape[2]; k++) {
2213 for (int l=0; l < shape[3]; l++) {
2214 value[i][j][k][l] = valueView(i,j,k,l);
2215 }
2216 }
2217 }
2218 }
2219 }
2220
2221 }
2222
2223 void
2224 Data::archiveData(const std::string fileName)
2225 {
2226 cout << "Archiving Data object to: " << fileName << endl;
2227
2228 //
2229 // Determine type of this Data object
2230 int dataType = -1;
2231
2232 if (isEmpty()) {
2233 dataType = 0;
2234 cout << "\tdataType: DataEmpty" << endl;
2235 }
2236 if (isConstant()) {
2237 dataType = 1;
2238 cout << "\tdataType: DataConstant" << endl;
2239 }
2240 if (isTagged()) {
2241 dataType = 2;
2242 cout << "\tdataType: DataTagged" << endl;
2243 }
2244 if (isExpanded()) {
2245 dataType = 3;
2246 cout << "\tdataType: DataExpanded" << endl;
2247 }
2248
2249 if (dataType == -1) {
2250 throw DataException("archiveData Error: undefined dataType");
2251 }
2252
2253 //
2254 // Collect data items common to all Data types
2255 int noSamples = getNumSamples();
2256 int noDPPSample = getNumDataPointsPerSample();
2257 int functionSpaceType = getFunctionSpace().getTypeCode();
2258 int dataPointRank = getDataPointRank();
2259 int dataPointSize = getDataPointSize();
2260 int dataLength = getLength();
2261 DataArrayView::ShapeType dataPointShape = getDataPointShape();
2262 vector<int> referenceNumbers(noSamples);
2263 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2264 referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
2265 }
2266 vector<int> tagNumbers(noSamples);
2267 if (isTagged()) {
2268 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2269 tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
2270 }
2271 }
2272
2273 cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2274 cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2275 cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2276
2277 //
2278 // Flatten Shape to an array of integers suitable for writing to file
2279 int flatShape[4] = {0,0,0,0};
2280 cout << "\tshape: < ";
2281 for (int dim=0; dim<dataPointRank; dim++) {
2282 flatShape[dim] = dataPointShape[dim];
2283 cout << dataPointShape[dim] << " ";
2284 }
2285 cout << ">" << endl;
2286
2287 //
2288 // Open archive file
2289 ofstream archiveFile;
2290 archiveFile.open(fileName.data(), ios::out);
2291
2292 if (!archiveFile.good()) {
2293 throw DataException("archiveData Error: problem opening archive file");
2294 }
2295
2296 //
2297 // Write common data items to archive file
2298 archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
2299 archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
2300 archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2301 archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2302 archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2303 archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2304 archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
2305 for (int dim = 0; dim < 4; dim++) {
2306 archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2307 }
2308 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2309 archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2310 }
2311 if (isTagged()) {
2312 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2313 archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2314 }
2315 }
2316
2317 if (!archiveFile.good()) {
2318 throw DataException("archiveData Error: problem writing to archive file");
2319 }
2320
2321 //
2322 // Archive underlying data values for each Data type
2323 int noValues;
2324 switch (dataType) {
2325 case 0:
2326 // DataEmpty
2327 noValues = 0;
2328 archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2329 cout << "\tnoValues: " << noValues << endl;
2330 break;
2331 case 1:
2332 // DataConstant
2333 noValues = m_data->getLength();
2334 archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2335 cout << "\tnoValues: " << noValues << endl;
2336 if (m_data->archiveData(archiveFile,noValues)) {
2337 throw DataException("archiveData Error: problem writing data to archive file");
2338 }
2339 break;
2340 case 2:
2341 // DataTagged
2342 noValues = m_data->getLength();
2343 archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2344 cout << "\tnoValues: " << noValues << endl;
2345 if (m_data->archiveData(archiveFile,noValues)) {
2346 throw DataException("archiveData Error: problem writing data to archive file");
2347 }
2348 break;
2349 case 3:
2350 // DataExpanded
2351 noValues = m_data->getLength();
2352 archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2353 cout << "\tnoValues: " << noValues << endl;
2354 if (m_data->archiveData(archiveFile,noValues)) {
2355 throw DataException("archiveData Error: problem writing data to archive file");
2356 }
2357 break;
2358 }
2359
2360 if (!archiveFile.good()) {
2361 throw DataException("archiveData Error: problem writing data to archive file");
2362 }
2363
2364 //
2365 // Close archive file
2366 archiveFile.close();
2367
2368 if (!archiveFile.good()) {
2369 throw DataException("archiveData Error: problem closing archive file");
2370 }
2371
2372 }
2373
2374 void
2375 Data::extractData(const std::string fileName,
2376 const FunctionSpace& fspace)
2377 {
2378 //
2379 // Can only extract Data to an object which is initially DataEmpty
2380 if (!isEmpty()) {
2381 throw DataException("extractData Error: can only extract to DataEmpty object");
2382 }
2383
2384 cout << "Extracting Data object from: " << fileName << endl;
2385
2386 int dataType;
2387 int noSamples;
2388 int noDPPSample;
2389 int functionSpaceType;
2390 int dataPointRank;
2391 int dataPointSize;
2392 int dataLength;
2393 DataArrayView::ShapeType dataPointShape;
2394 int flatShape[4];
2395
2396 //
2397 // Open the archive file
2398 ifstream archiveFile;
2399 archiveFile.open(fileName.data(), ios::in);
2400
2401 if (!archiveFile.good()) {
2402 throw DataException("extractData Error: problem opening archive file");
2403 }
2404
2405 //
2406 // Read common data items from archive file
2407 archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
2408 archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2409 archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2410 archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2411 archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2412 archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2413 archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
2414 for (int dim = 0; dim < 4; dim++) {
2415 archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2416 if (flatShape[dim]>0) {
2417 dataPointShape.push_back(flatShape[dim]);
2418 }
2419 }
2420 vector<int> referenceNumbers(noSamples);
2421 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2422 archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2423 }
2424 vector<int> tagNumbers(noSamples);
2425 if (dataType==2) {
2426 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2427 archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2428 }
2429 }
2430
2431 if (!archiveFile.good()) {
2432 throw DataException("extractData Error: problem reading from archive file");
2433 }
2434
2435 //
2436 // Verify the values just read from the archive file
2437 switch (dataType) {
2438 case 0:
2439 cout << "\tdataType: DataEmpty" << endl;
2440 break;
2441 case 1:
2442 cout << "\tdataType: DataConstant" << endl;
2443 break;
2444 case 2:
2445 cout << "\tdataType: DataTagged" << endl;
2446 break;
2447 case 3:
2448 cout << "\tdataType: DataExpanded" << endl;
2449 break;
2450 default:
2451 throw DataException("extractData Error: undefined dataType read from archive file");
2452 break;
2453 }
2454
2455 cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2456 cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2457 cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2458 cout << "\tshape: < ";
2459 for (int dim = 0; dim < dataPointRank; dim++) {
2460 cout << dataPointShape[dim] << " ";
2461 }
2462 cout << ">" << endl;
2463
2464 //
2465 // Verify that supplied FunctionSpace object is compatible with this Data object.
2466 if ( (fspace.getTypeCode()!=functionSpaceType) ||
2467 (fspace.getNumSamples()!=noSamples) ||
2468 (fspace.getNumDPPSample()!=noDPPSample)
2469 ) {
2470 throw DataException("extractData Error: incompatible FunctionSpace");
2471 }
2472 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2473 if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2474 throw DataException("extractData Error: incompatible FunctionSpace");
2475 }
2476 }
2477 if (dataType==2) {
2478 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2479 if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2480 throw DataException("extractData Error: incompatible FunctionSpace");
2481 }
2482 }
2483 }
2484
2485 //
2486 // Construct a DataVector to hold underlying data values
2487 DataVector dataVec(dataLength);
2488
2489 //
2490 // Load this DataVector with the appropriate values
2491 int noValues;
2492 archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2493 cout << "\tnoValues: " << noValues << endl;
2494 switch (dataType) {
2495 case 0:
2496 // DataEmpty
2497 if (noValues != 0) {
2498 throw DataException("extractData Error: problem reading data from archive file");
2499 }
2500 break;
2501 case 1:
2502 // DataConstant
2503 if (dataVec.extractData(archiveFile,noValues)) {
2504 throw DataException("extractData Error: problem reading data from archive file");
2505 }
2506 break;
2507 case 2:
2508 // DataTagged
2509 if (dataVec.extractData(archiveFile,noValues)) {
2510 throw DataException("extractData Error: problem reading data from archive file");
2511 }
2512 break;
2513 case 3:
2514 // DataExpanded
2515 if (dataVec.extractData(archiveFile,noValues)) {
2516 throw DataException("extractData Error: problem reading data from archive file");
2517 }
2518 break;
2519 }
2520
2521 if (!archiveFile.good()) {
2522 throw DataException("extractData Error: problem reading from archive file");
2523 }
2524
2525 //
2526 // Close archive file
2527 archiveFile.close();
2528
2529 if (!archiveFile.good()) {
2530 throw DataException("extractData Error: problem closing archive file");
2531 }
2532
2533 //
2534 // Construct an appropriate Data object
2535 DataAbstract* tempData;
2536 switch (dataType) {
2537 case 0:
2538 // DataEmpty
2539 tempData=new DataEmpty();
2540 break;
2541 case 1:
2542 // DataConstant
2543 tempData=new DataConstant(fspace,dataPointShape,dataVec);
2544 break;
2545 case 2:
2546 // DataTagged
2547 tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2548 break;
2549 case 3:
2550 // DataExpanded
2551 tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2552 break;
2553 }
2554 shared_ptr<DataAbstract> temp_data(tempData);
2555 m_data=temp_data;
2556 }
2557
2558 ostream& escript::operator<<(ostream& o, const Data& data)
2559 {
2560 o << data.toString();
2561 return o;
2562 }
2563
2564 /* Member functions specific to the MPI implementation */
2565
2566 void
2567 Data::print()
2568 {
2569 int i,j;
2570
2571 printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2572 for( i=0; i<getNumSamples(); i++ )
2573 {
2574 printf( "[%6d]", i );
2575 for( j=0; j<getNumDataPointsPerSample(); j++ )
2576 printf( "\t%10.7g", (getSampleData(i))[j] );
2577 printf( "\n" );
2578 }
2579 }
2580
2581 int
2582 Data::get_MPISize() const
2583 {
2584 int error, size;
2585 #ifdef PASO_MPI
2586 error = MPI_Comm_size( get_MPIComm(), &size );
2587 #else
2588 size = 1;
2589 #endif
2590 return size;
2591 }
2592
2593 int
2594 Data::get_MPIRank() const
2595 {
2596 int error, rank;
2597 #ifdef PASO_MPI
2598 error = MPI_Comm_rank( get_MPIComm(), &rank );
2599 #else
2600 rank = 0;
2601 #endif
2602 return rank;
2603 }
2604
2605 MPI_Comm
2606 Data::get_MPIComm() const
2607 {
2608 #ifdef PASO_MPI
2609 return MPI_COMM_WORLD;
2610 #else
2611 return -1;
2612 #endif
2613 }
2614

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26