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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 790 - (show annotations)
Wed Jul 26 23:12:34 2006 UTC (12 years, 10 months ago) by bcumming
File size: 64492 byte(s)
changes to escript/py_src/pdetools.py and /escript/src/Data.h/.cpp to
make the Locator work in MPI. escript::Data::mindp now returns a 3 tuple,
with the MPI rank of the process on which the minimum value occurs
included. escript::Data::convertToNumArrayFromDPNo also takes the ProcNo
to perform the MPI reduction.

This had to be implemented in both the MPI and non-MPI versions to allow
the necesary changes to the Python code in pdetools.py. In the non-MPI
version ProcNo is set to 0. This works for the explicit scripts tested
thus far, however if it causes problems in your scripts contact Ben or
Lutz, or revert the three files (pdetools.py, Data.h and Data.cpp) to
the previous version.  


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::trace() const
1320 {
1321 #if defined DOPROF
1322 profData->reduction2++;
1323 #endif
1324 Trace trace_func;
1325 return dp_algorithm(trace_func,0);
1326 }
1327
1328 Data
1329 Data::symmetric() const
1330 {
1331 #if defined DOPROF
1332 profData->unary++;
1333 #endif
1334 // check input
1335 DataArrayView::ShapeType s=getDataPointShape();
1336 if (getDataPointRank()==2) {
1337 if(s[0] != s[1])
1338 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1339 }
1340 else if (getDataPointRank()==4) {
1341 if(!(s[0] == s[2] && s[1] == s[3]))
1342 throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1343 }
1344 else {
1345 throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1346 }
1347 Data ev(0.,getDataPointShape(),getFunctionSpace());
1348 ev.typeMatchRight(*this);
1349 m_data->symmetric(ev.m_data.get());
1350 return ev;
1351 }
1352
1353 Data
1354 Data::nonsymmetric() const
1355 {
1356 #if defined DOPROF
1357 profData->unary++;
1358 #endif
1359 // check input
1360 DataArrayView::ShapeType s=getDataPointShape();
1361 if (getDataPointRank()==2) {
1362 if(s[0] != s[1])
1363 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1364 DataArrayView::ShapeType ev_shape;
1365 ev_shape.push_back(s[0]);
1366 ev_shape.push_back(s[1]);
1367 Data ev(0.,ev_shape,getFunctionSpace());
1368 ev.typeMatchRight(*this);
1369 m_data->nonsymmetric(ev.m_data.get());
1370 return ev;
1371 }
1372 else if (getDataPointRank()==4) {
1373 if(!(s[0] == s[2] && s[1] == s[3]))
1374 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1375 DataArrayView::ShapeType ev_shape;
1376 ev_shape.push_back(s[0]);
1377 ev_shape.push_back(s[1]);
1378 ev_shape.push_back(s[2]);
1379 ev_shape.push_back(s[3]);
1380 Data ev(0.,ev_shape,getFunctionSpace());
1381 ev.typeMatchRight(*this);
1382 m_data->nonsymmetric(ev.m_data.get());
1383 return ev;
1384 }
1385 else {
1386 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1387 }
1388 }
1389
1390 Data
1391 Data::matrixtrace(int axis_offset) const
1392 {
1393 #if defined DOPROF
1394 profData->unary++;
1395 #endif
1396 DataArrayView::ShapeType s=getDataPointShape();
1397 if (getDataPointRank()==2) {
1398 DataArrayView::ShapeType ev_shape;
1399 Data ev(0.,ev_shape,getFunctionSpace());
1400 ev.typeMatchRight(*this);
1401 m_data->matrixtrace(ev.m_data.get(), axis_offset);
1402 return ev;
1403 }
1404 if (getDataPointRank()==3) {
1405 DataArrayView::ShapeType ev_shape;
1406 if (axis_offset==0) {
1407 int s2=s[2];
1408 ev_shape.push_back(s2);
1409 }
1410 else if (axis_offset==1) {
1411 int s0=s[0];
1412 ev_shape.push_back(s0);
1413 }
1414 Data ev(0.,ev_shape,getFunctionSpace());
1415 ev.typeMatchRight(*this);
1416 m_data->matrixtrace(ev.m_data.get(), axis_offset);
1417 return ev;
1418 }
1419 if (getDataPointRank()==4) {
1420 DataArrayView::ShapeType ev_shape;
1421 if (axis_offset==0) {
1422 ev_shape.push_back(s[2]);
1423 ev_shape.push_back(s[3]);
1424 }
1425 else if (axis_offset==1) {
1426 ev_shape.push_back(s[0]);
1427 ev_shape.push_back(s[3]);
1428 }
1429 else if (axis_offset==2) {
1430 ev_shape.push_back(s[0]);
1431 ev_shape.push_back(s[1]);
1432 }
1433 Data ev(0.,ev_shape,getFunctionSpace());
1434 ev.typeMatchRight(*this);
1435 m_data->matrixtrace(ev.m_data.get(), axis_offset);
1436 return ev;
1437 }
1438 else {
1439 throw DataException("Error - Data::matrixtrace can only be calculated for rank 2, 3 or 4 object.");
1440 }
1441 }
1442
1443 Data
1444 Data::transpose(int axis_offset) const
1445 {
1446 #if defined DOPROF
1447 profData->reduction2++;
1448 #endif
1449 DataArrayView::ShapeType s=getDataPointShape();
1450 DataArrayView::ShapeType ev_shape;
1451 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1452 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1453 int rank=getDataPointRank();
1454 if (axis_offset<0 || axis_offset>rank) {
1455 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1456 }
1457 for (int i=0; i<rank; i++) {
1458 int index = (axis_offset+i)%rank;
1459 ev_shape.push_back(s[index]); // Append to new shape
1460 }
1461 Data ev(0.,ev_shape,getFunctionSpace());
1462 ev.typeMatchRight(*this);
1463 m_data->transpose(ev.m_data.get(), axis_offset);
1464 return ev;
1465 }
1466
1467 Data
1468 Data::eigenvalues() const
1469 {
1470 #if defined DOPROF
1471 profData->unary++;
1472 #endif
1473 // check input
1474 DataArrayView::ShapeType s=getDataPointShape();
1475 if (getDataPointRank()!=2)
1476 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1477 if(s[0] != s[1])
1478 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1479 // create return
1480 DataArrayView::ShapeType ev_shape(1,s[0]);
1481 Data ev(0.,ev_shape,getFunctionSpace());
1482 ev.typeMatchRight(*this);
1483 m_data->eigenvalues(ev.m_data.get());
1484 return ev;
1485 }
1486
1487 const boost::python::tuple
1488 Data::eigenvalues_and_eigenvectors(const double tol) const
1489 {
1490 #if defined DOPROF
1491 profData->unary++;
1492 #endif
1493 DataArrayView::ShapeType s=getDataPointShape();
1494 if (getDataPointRank()!=2)
1495 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1496 if(s[0] != s[1])
1497 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1498 // create return
1499 DataArrayView::ShapeType ev_shape(1,s[0]);
1500 Data ev(0.,ev_shape,getFunctionSpace());
1501 ev.typeMatchRight(*this);
1502 DataArrayView::ShapeType V_shape(2,s[0]);
1503 Data V(0.,V_shape,getFunctionSpace());
1504 V.typeMatchRight(*this);
1505 m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1506 return make_tuple(boost::python::object(ev),boost::python::object(V));
1507 }
1508
1509 const boost::python::tuple
1510 Data::mindp() const
1511 {
1512 // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1513 // abort (for unknown reasons) if there are openmp directives with it in the
1514 // surrounding function
1515
1516 int SampleNo;
1517 int DataPointNo;
1518 int ProcNo;
1519 calc_mindp(ProcNo,SampleNo,DataPointNo);
1520 return make_tuple(ProcNo,SampleNo,DataPointNo);
1521 }
1522
1523 void
1524 Data::calc_mindp( int& ProcNo,
1525 int& SampleNo,
1526 int& DataPointNo) const
1527 {
1528 int i,j;
1529 int lowi=0,lowj=0;
1530 double min=numeric_limits<double>::max();
1531
1532 Data temp=minval();
1533
1534 int numSamples=temp.getNumSamples();
1535 int numDPPSample=temp.getNumDataPointsPerSample();
1536
1537 double next,local_min;
1538 int local_lowi,local_lowj;
1539
1540 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1541 {
1542 local_min=min;
1543 #pragma omp for private(i,j) schedule(static)
1544 for (i=0; i<numSamples; i++) {
1545 for (j=0; j<numDPPSample; j++) {
1546 next=temp.getDataPoint(i,j)();
1547 if (next<local_min) {
1548 local_min=next;
1549 local_lowi=i;
1550 local_lowj=j;
1551 }
1552 }
1553 }
1554 #pragma omp critical
1555 if (local_min<min) {
1556 min=local_min;
1557 lowi=local_lowi;
1558 lowj=local_lowj;
1559 }
1560 }
1561
1562 #ifdef PASO_MPI
1563 // determine the processor on which the minimum occurs
1564 next = temp.getDataPoint(lowi,lowj)();
1565 int lowProc = 0;
1566 double *globalMins = new double[get_MPISize()+1];
1567 int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1568
1569 if( get_MPIRank()==0 ){
1570 next = globalMins[lowProc];
1571 for( i=1; i<get_MPISize(); i++ )
1572 if( next>globalMins[i] ){
1573 lowProc = i;
1574 next = globalMins[i];
1575 }
1576 }
1577 MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1578
1579 delete [] globalMins;
1580 ProcNo = lowProc;
1581 #else
1582 ProcNo = 0;
1583 #endif
1584 SampleNo = lowi;
1585 DataPointNo = lowj;
1586 }
1587
1588 void
1589 Data::saveDX(std::string fileName) const
1590 {
1591 boost::python::dict args;
1592 args["data"]=boost::python::object(this);
1593 getDomain().saveDX(fileName,args);
1594 return;
1595 }
1596
1597 void
1598 Data::saveVTK(std::string fileName) const
1599 {
1600 boost::python::dict args;
1601 args["data"]=boost::python::object(this);
1602 getDomain().saveVTK(fileName,args);
1603 return;
1604 }
1605
1606 Data&
1607 Data::operator+=(const Data& right)
1608 {
1609 if (isProtected()) {
1610 throw DataException("Error - attempt to update protected Data object.");
1611 }
1612 #if defined DOPROF
1613 profData->binary++;
1614 #endif
1615 binaryOp(right,plus<double>());
1616 return (*this);
1617 }
1618
1619 Data&
1620 Data::operator+=(const boost::python::object& right)
1621 {
1622 if (isProtected()) {
1623 throw DataException("Error - attempt to update protected Data object.");
1624 }
1625 #if defined DOPROF
1626 profData->binary++;
1627 #endif
1628 binaryOp(right,plus<double>());
1629 return (*this);
1630 }
1631
1632 Data&
1633 Data::operator-=(const Data& right)
1634 {
1635 if (isProtected()) {
1636 throw DataException("Error - attempt to update protected Data object.");
1637 }
1638 #if defined DOPROF
1639 profData->binary++;
1640 #endif
1641 binaryOp(right,minus<double>());
1642 return (*this);
1643 }
1644
1645 Data&
1646 Data::operator-=(const boost::python::object& right)
1647 {
1648 if (isProtected()) {
1649 throw DataException("Error - attempt to update protected Data object.");
1650 }
1651 #if defined DOPROF
1652 profData->binary++;
1653 #endif
1654 binaryOp(right,minus<double>());
1655 return (*this);
1656 }
1657
1658 Data&
1659 Data::operator*=(const Data& right)
1660 {
1661 if (isProtected()) {
1662 throw DataException("Error - attempt to update protected Data object.");
1663 }
1664 #if defined DOPROF
1665 profData->binary++;
1666 #endif
1667 binaryOp(right,multiplies<double>());
1668 return (*this);
1669 }
1670
1671 Data&
1672 Data::operator*=(const boost::python::object& right)
1673 {
1674 if (isProtected()) {
1675 throw DataException("Error - attempt to update protected Data object.");
1676 }
1677 #if defined DOPROF
1678 profData->binary++;
1679 #endif
1680 binaryOp(right,multiplies<double>());
1681 return (*this);
1682 }
1683
1684 Data&
1685 Data::operator/=(const Data& right)
1686 {
1687 if (isProtected()) {
1688 throw DataException("Error - attempt to update protected Data object.");
1689 }
1690 #if defined DOPROF
1691 profData->binary++;
1692 #endif
1693 binaryOp(right,divides<double>());
1694 return (*this);
1695 }
1696
1697 Data&
1698 Data::operator/=(const boost::python::object& right)
1699 {
1700 if (isProtected()) {
1701 throw DataException("Error - attempt to update protected Data object.");
1702 }
1703 #if defined DOPROF
1704 profData->binary++;
1705 #endif
1706 binaryOp(right,divides<double>());
1707 return (*this);
1708 }
1709
1710 Data
1711 Data::rpowO(const boost::python::object& left) const
1712 {
1713 if (isProtected()) {
1714 throw DataException("Error - attempt to update protected Data object.");
1715 }
1716 #if defined DOPROF
1717 profData->binary++;
1718 #endif
1719 Data left_d(left,*this);
1720 return left_d.powD(*this);
1721 }
1722
1723 Data
1724 Data::powO(const boost::python::object& right) const
1725 {
1726 #if defined DOPROF
1727 profData->binary++;
1728 #endif
1729 Data result;
1730 result.copy(*this);
1731 result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1732 return result;
1733 }
1734
1735 Data
1736 Data::powD(const Data& right) const
1737 {
1738 #if defined DOPROF
1739 profData->binary++;
1740 #endif
1741 Data result;
1742 result.copy(*this);
1743 result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1744 return result;
1745 }
1746
1747
1748 //
1749 // NOTE: It is essential to specify the namespace this operator belongs to
1750 Data
1751 escript::operator+(const Data& left, const Data& right)
1752 {
1753 Data result;
1754 //
1755 // perform a deep copy
1756 result.copy(left);
1757 result+=right;
1758 return result;
1759 }
1760
1761 //
1762 // NOTE: It is essential to specify the namespace this operator belongs to
1763 Data
1764 escript::operator-(const Data& left, const Data& right)
1765 {
1766 Data result;
1767 //
1768 // perform a deep copy
1769 result.copy(left);
1770 result-=right;
1771 return result;
1772 }
1773
1774 //
1775 // NOTE: It is essential to specify the namespace this operator belongs to
1776 Data
1777 escript::operator*(const Data& left, const Data& right)
1778 {
1779 Data result;
1780 //
1781 // perform a deep copy
1782 result.copy(left);
1783 result*=right;
1784 return result;
1785 }
1786
1787 //
1788 // NOTE: It is essential to specify the namespace this operator belongs to
1789 Data
1790 escript::operator/(const Data& left, const Data& right)
1791 {
1792 Data result;
1793 //
1794 // perform a deep copy
1795 result.copy(left);
1796 result/=right;
1797 return result;
1798 }
1799
1800 //
1801 // NOTE: It is essential to specify the namespace this operator belongs to
1802 Data
1803 escript::operator+(const Data& left, const boost::python::object& right)
1804 {
1805 //
1806 // Convert to DataArray format if possible
1807 DataArray temp(right);
1808 Data result;
1809 //
1810 // perform a deep copy
1811 result.copy(left);
1812 result+=right;
1813 return result;
1814 }
1815
1816 //
1817 // NOTE: It is essential to specify the namespace this operator belongs to
1818 Data
1819 escript::operator-(const Data& left, const boost::python::object& right)
1820 {
1821 //
1822 // Convert to DataArray format if possible
1823 DataArray temp(right);
1824 Data result;
1825 //
1826 // perform a deep copy
1827 result.copy(left);
1828 result-=right;
1829 return result;
1830 }
1831
1832 //
1833 // NOTE: It is essential to specify the namespace this operator belongs to
1834 Data
1835 escript::operator*(const Data& left, const boost::python::object& right)
1836 {
1837 //
1838 // Convert to DataArray format if possible
1839 DataArray temp(right);
1840 Data result;
1841 //
1842 // perform a deep copy
1843 result.copy(left);
1844 result*=right;
1845 return result;
1846 }
1847
1848 //
1849 // NOTE: It is essential to specify the namespace this operator belongs to
1850 Data
1851 escript::operator/(const Data& left, const boost::python::object& right)
1852 {
1853 //
1854 // Convert to DataArray format if possible
1855 DataArray temp(right);
1856 Data result;
1857 //
1858 // perform a deep copy
1859 result.copy(left);
1860 result/=right;
1861 return result;
1862 }
1863
1864 //
1865 // NOTE: It is essential to specify the namespace this operator belongs to
1866 Data
1867 escript::operator+(const boost::python::object& left, const Data& right)
1868 {
1869 //
1870 // Construct the result using the given value and the other parameters
1871 // from right
1872 Data result(left,right);
1873 result+=right;
1874 return result;
1875 }
1876
1877 //
1878 // NOTE: It is essential to specify the namespace this operator belongs to
1879 Data
1880 escript::operator-(const boost::python::object& left, const Data& right)
1881 {
1882 //
1883 // Construct the result using the given value and the other parameters
1884 // from right
1885 Data result(left,right);
1886 result-=right;
1887 return result;
1888 }
1889
1890 //
1891 // NOTE: It is essential to specify the namespace this operator belongs to
1892 Data
1893 escript::operator*(const boost::python::object& left, const Data& right)
1894 {
1895 //
1896 // Construct the result using the given value and the other parameters
1897 // from right
1898 Data result(left,right);
1899 result*=right;
1900 return result;
1901 }
1902
1903 //
1904 // NOTE: It is essential to specify the namespace this operator belongs to
1905 Data
1906 escript::operator/(const boost::python::object& left, const Data& right)
1907 {
1908 //
1909 // Construct the result using the given value and the other parameters
1910 // from right
1911 Data result(left,right);
1912 result/=right;
1913 return result;
1914 }
1915
1916 //
1917 //bool escript::operator==(const Data& left, const Data& right)
1918 //{
1919 // /*
1920 // NB: this operator does very little at this point, and isn't to
1921 // be relied on. Requires further implementation.
1922 // */
1923 //
1924 // bool ret;
1925 //
1926 // if (left.isEmpty()) {
1927 // if(!right.isEmpty()) {
1928 // ret = false;
1929 // } else {
1930 // ret = true;
1931 // }
1932 // }
1933 //
1934 // if (left.isConstant()) {
1935 // if(!right.isConstant()) {
1936 // ret = false;
1937 // } else {
1938 // ret = true;
1939 // }
1940 // }
1941 //
1942 // if (left.isTagged()) {
1943 // if(!right.isTagged()) {
1944 // ret = false;
1945 // } else {
1946 // ret = true;
1947 // }
1948 // }
1949 //
1950 // if (left.isExpanded()) {
1951 // if(!right.isExpanded()) {
1952 // ret = false;
1953 // } else {
1954 // ret = true;
1955 // }
1956 // }
1957 //
1958 // return ret;
1959 //}
1960
1961 /* TODO */
1962 /* global reduction */
1963 Data
1964 Data::getItem(const boost::python::object& key) const
1965 {
1966 const DataArrayView& view=getPointDataView();
1967
1968 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1969
1970 if (slice_region.size()!=view.getRank()) {
1971 throw DataException("Error - slice size does not match Data rank.");
1972 }
1973
1974 return getSlice(slice_region);
1975 }
1976
1977 /* TODO */
1978 /* global reduction */
1979 Data
1980 Data::getSlice(const DataArrayView::RegionType& region) const
1981 {
1982 #if defined DOPROF
1983 profData->slicing++;
1984 #endif
1985 return Data(*this,region);
1986 }
1987
1988 /* TODO */
1989 /* global reduction */
1990 void
1991 Data::setItemO(const boost::python::object& key,
1992 const boost::python::object& value)
1993 {
1994 Data tempData(value,getFunctionSpace());
1995 setItemD(key,tempData);
1996 }
1997
1998 /* TODO */
1999 /* global reduction */
2000 void
2001 Data::setItemD(const boost::python::object& key,
2002 const Data& value)
2003 {
2004 const DataArrayView& view=getPointDataView();
2005
2006 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
2007 if (slice_region.size()!=view.getRank()) {
2008 throw DataException("Error - slice size does not match Data rank.");
2009 }
2010 if (getFunctionSpace()!=value.getFunctionSpace()) {
2011 setSlice(Data(value,getFunctionSpace()),slice_region);
2012 } else {
2013 setSlice(value,slice_region);
2014 }
2015 }
2016
2017 /* TODO */
2018 /* global reduction */
2019 void
2020 Data::setSlice(const Data& value,
2021 const DataArrayView::RegionType& region)
2022 {
2023 if (isProtected()) {
2024 throw DataException("Error - attempt to update protected Data object.");
2025 }
2026 #if defined DOPROF
2027 profData->slicing++;
2028 #endif
2029 Data tempValue(value);
2030 typeMatchLeft(tempValue);
2031 typeMatchRight(tempValue);
2032 m_data->setSlice(tempValue.m_data.get(),region);
2033 }
2034
2035 void
2036 Data::typeMatchLeft(Data& right) const
2037 {
2038 if (isExpanded()){
2039 right.expand();
2040 } else if (isTagged()) {
2041 if (right.isConstant()) {
2042 right.tag();
2043 }
2044 }
2045 }
2046
2047 void
2048 Data::typeMatchRight(const Data& right)
2049 {
2050 if (isTagged()) {
2051 if (right.isExpanded()) {
2052 expand();
2053 }
2054 } else if (isConstant()) {
2055 if (right.isExpanded()) {
2056 expand();
2057 } else if (right.isTagged()) {
2058 tag();
2059 }
2060 }
2061 }
2062
2063 /* TODO */
2064 /* global reduction */
2065 void
2066 Data::setTaggedValue(int tagKey,
2067 const boost::python::object& value)
2068 {
2069 if (isProtected()) {
2070 throw DataException("Error - attempt to update protected Data object.");
2071 }
2072 //
2073 // Ensure underlying data object is of type DataTagged
2074 tag();
2075
2076 if (!isTagged()) {
2077 throw DataException("Error - DataTagged conversion failed!!");
2078 }
2079
2080 //
2081 // Construct DataArray from boost::python::object input value
2082 DataArray valueDataArray(value);
2083
2084 //
2085 // Call DataAbstract::setTaggedValue
2086 m_data->setTaggedValue(tagKey,valueDataArray.getView());
2087 }
2088
2089 /* TODO */
2090 /* global reduction */
2091 void
2092 Data::setTaggedValueFromCPP(int tagKey,
2093 const DataArrayView& value)
2094 {
2095 if (isProtected()) {
2096 throw DataException("Error - attempt to update protected Data object.");
2097 }
2098 //
2099 // Ensure underlying data object is of type DataTagged
2100 tag();
2101
2102 if (!isTagged()) {
2103 throw DataException("Error - DataTagged conversion failed!!");
2104 }
2105
2106 //
2107 // Call DataAbstract::setTaggedValue
2108 m_data->setTaggedValue(tagKey,value);
2109 }
2110
2111 /* TODO */
2112 /* global reduction */
2113 int
2114 Data::getTagNumber(int dpno)
2115 {
2116 return m_data->getTagNumber(dpno);
2117 }
2118
2119 /* TODO */
2120 /* global reduction */
2121 void
2122 Data::setRefValue(int ref,
2123 const boost::python::numeric::array& value)
2124 {
2125 if (isProtected()) {
2126 throw DataException("Error - attempt to update protected Data object.");
2127 }
2128 //
2129 // Construct DataArray from boost::python::object input value
2130 DataArray valueDataArray(value);
2131
2132 //
2133 // Call DataAbstract::setRefValue
2134 m_data->setRefValue(ref,valueDataArray);
2135 }
2136
2137 /* TODO */
2138 /* global reduction */
2139 void
2140 Data::getRefValue(int ref,
2141 boost::python::numeric::array& value)
2142 {
2143 //
2144 // Construct DataArray for boost::python::object return value
2145 DataArray valueDataArray(value);
2146
2147 //
2148 // Load DataArray with values from data-points specified by ref
2149 m_data->getRefValue(ref,valueDataArray);
2150
2151 //
2152 // Load values from valueDataArray into return numarray
2153
2154 // extract the shape of the numarray
2155 int rank = value.getrank();
2156 DataArrayView::ShapeType shape;
2157 for (int i=0; i < rank; i++) {
2158 shape.push_back(extract<int>(value.getshape()[i]));
2159 }
2160
2161 // and load the numarray with the data from the DataArray
2162 DataArrayView valueView = valueDataArray.getView();
2163
2164 if (rank==0) {
2165 boost::python::numeric::array temp_numArray(valueView());
2166 value = temp_numArray;
2167 }
2168 if (rank==1) {
2169 for (int i=0; i < shape[0]; i++) {
2170 value[i] = valueView(i);
2171 }
2172 }
2173 if (rank==2) {
2174 for (int i=0; i < shape[0]; i++) {
2175 for (int j=0; j < shape[1]; j++) {
2176 value[i][j] = valueView(i,j);
2177 }
2178 }
2179 }
2180 if (rank==3) {
2181 for (int i=0; i < shape[0]; i++) {
2182 for (int j=0; j < shape[1]; j++) {
2183 for (int k=0; k < shape[2]; k++) {
2184 value[i][j][k] = valueView(i,j,k);
2185 }
2186 }
2187 }
2188 }
2189 if (rank==4) {
2190 for (int i=0; i < shape[0]; i++) {
2191 for (int j=0; j < shape[1]; j++) {
2192 for (int k=0; k < shape[2]; k++) {
2193 for (int l=0; l < shape[3]; l++) {
2194 value[i][j][k][l] = valueView(i,j,k,l);
2195 }
2196 }
2197 }
2198 }
2199 }
2200
2201 }
2202
2203 void
2204 Data::archiveData(const std::string fileName)
2205 {
2206 cout << "Archiving Data object to: " << fileName << endl;
2207
2208 //
2209 // Determine type of this Data object
2210 int dataType = -1;
2211
2212 if (isEmpty()) {
2213 dataType = 0;
2214 cout << "\tdataType: DataEmpty" << endl;
2215 }
2216 if (isConstant()) {
2217 dataType = 1;
2218 cout << "\tdataType: DataConstant" << endl;
2219 }
2220 if (isTagged()) {
2221 dataType = 2;
2222 cout << "\tdataType: DataTagged" << endl;
2223 }
2224 if (isExpanded()) {
2225 dataType = 3;
2226 cout << "\tdataType: DataExpanded" << endl;
2227 }
2228
2229 if (dataType == -1) {
2230 throw DataException("archiveData Error: undefined dataType");
2231 }
2232
2233 //
2234 // Collect data items common to all Data types
2235 int noSamples = getNumSamples();
2236 int noDPPSample = getNumDataPointsPerSample();
2237 int functionSpaceType = getFunctionSpace().getTypeCode();
2238 int dataPointRank = getDataPointRank();
2239 int dataPointSize = getDataPointSize();
2240 int dataLength = getLength();
2241 DataArrayView::ShapeType dataPointShape = getDataPointShape();
2242 vector<int> referenceNumbers(noSamples);
2243 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2244 referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
2245 }
2246 vector<int> tagNumbers(noSamples);
2247 if (isTagged()) {
2248 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2249 tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
2250 }
2251 }
2252
2253 cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2254 cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2255 cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2256
2257 //
2258 // Flatten Shape to an array of integers suitable for writing to file
2259 int flatShape[4] = {0,0,0,0};
2260 cout << "\tshape: < ";
2261 for (int dim=0; dim<dataPointRank; dim++) {
2262 flatShape[dim] = dataPointShape[dim];
2263 cout << dataPointShape[dim] << " ";
2264 }
2265 cout << ">" << endl;
2266
2267 //
2268 // Open archive file
2269 ofstream archiveFile;
2270 archiveFile.open(fileName.data(), ios::out);
2271
2272 if (!archiveFile.good()) {
2273 throw DataException("archiveData Error: problem opening archive file");
2274 }
2275
2276 //
2277 // Write common data items to archive file
2278 archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
2279 archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
2280 archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2281 archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2282 archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2283 archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2284 archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
2285 for (int dim = 0; dim < 4; dim++) {
2286 archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2287 }
2288 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2289 archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2290 }
2291 if (isTagged()) {
2292 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2293 archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2294 }
2295 }
2296
2297 if (!archiveFile.good()) {
2298 throw DataException("archiveData Error: problem writing to archive file");
2299 }
2300
2301 //
2302 // Archive underlying data values for each Data type
2303 int noValues;
2304 switch (dataType) {
2305 case 0:
2306 // DataEmpty
2307 noValues = 0;
2308 archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2309 cout << "\tnoValues: " << noValues << endl;
2310 break;
2311 case 1:
2312 // DataConstant
2313 noValues = m_data->getLength();
2314 archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2315 cout << "\tnoValues: " << noValues << endl;
2316 if (m_data->archiveData(archiveFile,noValues)) {
2317 throw DataException("archiveData Error: problem writing data to archive file");
2318 }
2319 break;
2320 case 2:
2321 // DataTagged
2322 noValues = m_data->getLength();
2323 archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2324 cout << "\tnoValues: " << noValues << endl;
2325 if (m_data->archiveData(archiveFile,noValues)) {
2326 throw DataException("archiveData Error: problem writing data to archive file");
2327 }
2328 break;
2329 case 3:
2330 // DataExpanded
2331 noValues = m_data->getLength();
2332 archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2333 cout << "\tnoValues: " << noValues << endl;
2334 if (m_data->archiveData(archiveFile,noValues)) {
2335 throw DataException("archiveData Error: problem writing data to archive file");
2336 }
2337 break;
2338 }
2339
2340 if (!archiveFile.good()) {
2341 throw DataException("archiveData Error: problem writing data to archive file");
2342 }
2343
2344 //
2345 // Close archive file
2346 archiveFile.close();
2347
2348 if (!archiveFile.good()) {
2349 throw DataException("archiveData Error: problem closing archive file");
2350 }
2351
2352 }
2353
2354 void
2355 Data::extractData(const std::string fileName,
2356 const FunctionSpace& fspace)
2357 {
2358 //
2359 // Can only extract Data to an object which is initially DataEmpty
2360 if (!isEmpty()) {
2361 throw DataException("extractData Error: can only extract to DataEmpty object");
2362 }
2363
2364 cout << "Extracting Data object from: " << fileName << endl;
2365
2366 int dataType;
2367 int noSamples;
2368 int noDPPSample;
2369 int functionSpaceType;
2370 int dataPointRank;
2371 int dataPointSize;
2372 int dataLength;
2373 DataArrayView::ShapeType dataPointShape;
2374 int flatShape[4];
2375
2376 //
2377 // Open the archive file
2378 ifstream archiveFile;
2379 archiveFile.open(fileName.data(), ios::in);
2380
2381 if (!archiveFile.good()) {
2382 throw DataException("extractData Error: problem opening archive file");
2383 }
2384
2385 //
2386 // Read common data items from archive file
2387 archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
2388 archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2389 archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2390 archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2391 archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2392 archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2393 archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
2394 for (int dim = 0; dim < 4; dim++) {
2395 archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2396 if (flatShape[dim]>0) {
2397 dataPointShape.push_back(flatShape[dim]);
2398 }
2399 }
2400 vector<int> referenceNumbers(noSamples);
2401 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2402 archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2403 }
2404 vector<int> tagNumbers(noSamples);
2405 if (dataType==2) {
2406 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2407 archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2408 }
2409 }
2410
2411 if (!archiveFile.good()) {
2412 throw DataException("extractData Error: problem reading from archive file");
2413 }
2414
2415 //
2416 // Verify the values just read from the archive file
2417 switch (dataType) {
2418 case 0:
2419 cout << "\tdataType: DataEmpty" << endl;
2420 break;
2421 case 1:
2422 cout << "\tdataType: DataConstant" << endl;
2423 break;
2424 case 2:
2425 cout << "\tdataType: DataTagged" << endl;
2426 break;
2427 case 3:
2428 cout << "\tdataType: DataExpanded" << endl;
2429 break;
2430 default:
2431 throw DataException("extractData Error: undefined dataType read from archive file");
2432 break;
2433 }
2434
2435 cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2436 cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2437 cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2438 cout << "\tshape: < ";
2439 for (int dim = 0; dim < dataPointRank; dim++) {
2440 cout << dataPointShape[dim] << " ";
2441 }
2442 cout << ">" << endl;
2443
2444 //
2445 // Verify that supplied FunctionSpace object is compatible with this Data object.
2446 if ( (fspace.getTypeCode()!=functionSpaceType) ||
2447 (fspace.getNumSamples()!=noSamples) ||
2448 (fspace.getNumDPPSample()!=noDPPSample)
2449 ) {
2450 throw DataException("extractData Error: incompatible FunctionSpace");
2451 }
2452 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2453 if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2454 throw DataException("extractData Error: incompatible FunctionSpace");
2455 }
2456 }
2457 if (dataType==2) {
2458 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2459 if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2460 throw DataException("extractData Error: incompatible FunctionSpace");
2461 }
2462 }
2463 }
2464
2465 //
2466 // Construct a DataVector to hold underlying data values
2467 DataVector dataVec(dataLength);
2468
2469 //
2470 // Load this DataVector with the appropriate values
2471 int noValues;
2472 archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2473 cout << "\tnoValues: " << noValues << endl;
2474 switch (dataType) {
2475 case 0:
2476 // DataEmpty
2477 if (noValues != 0) {
2478 throw DataException("extractData Error: problem reading data from archive file");
2479 }
2480 break;
2481 case 1:
2482 // DataConstant
2483 if (dataVec.extractData(archiveFile,noValues)) {
2484 throw DataException("extractData Error: problem reading data from archive file");
2485 }
2486 break;
2487 case 2:
2488 // DataTagged
2489 if (dataVec.extractData(archiveFile,noValues)) {
2490 throw DataException("extractData Error: problem reading data from archive file");
2491 }
2492 break;
2493 case 3:
2494 // DataExpanded
2495 if (dataVec.extractData(archiveFile,noValues)) {
2496 throw DataException("extractData Error: problem reading data from archive file");
2497 }
2498 break;
2499 }
2500
2501 if (!archiveFile.good()) {
2502 throw DataException("extractData Error: problem reading from archive file");
2503 }
2504
2505 //
2506 // Close archive file
2507 archiveFile.close();
2508
2509 if (!archiveFile.good()) {
2510 throw DataException("extractData Error: problem closing archive file");
2511 }
2512
2513 //
2514 // Construct an appropriate Data object
2515 DataAbstract* tempData;
2516 switch (dataType) {
2517 case 0:
2518 // DataEmpty
2519 tempData=new DataEmpty();
2520 break;
2521 case 1:
2522 // DataConstant
2523 tempData=new DataConstant(fspace,dataPointShape,dataVec);
2524 break;
2525 case 2:
2526 // DataTagged
2527 tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2528 break;
2529 case 3:
2530 // DataExpanded
2531 tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2532 break;
2533 }
2534 shared_ptr<DataAbstract> temp_data(tempData);
2535 m_data=temp_data;
2536 }
2537
2538 ostream& escript::operator<<(ostream& o, const Data& data)
2539 {
2540 o << data.toString();
2541 return o;
2542 }
2543
2544 /* Member functions specific to the MPI implementation */
2545
2546 void
2547 Data::print()
2548 {
2549 int i,j;
2550
2551 printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2552 for( i=0; i<getNumSamples(); i++ )
2553 {
2554 printf( "[%6d]", i );
2555 for( j=0; j<getNumDataPointsPerSample(); j++ )
2556 printf( "\t%10.7g", (getSampleData(i))[j] );
2557 printf( "\n" );
2558 }
2559 }
2560
2561 int
2562 Data::get_MPISize() const
2563 {
2564 int error, size;
2565 #ifdef PASO_MPI
2566 error = MPI_Comm_size( get_MPIComm(), &size );
2567 #else
2568 size = 1;
2569 #endif
2570 return size;
2571 }
2572
2573 int
2574 Data::get_MPIRank() const
2575 {
2576 int error, rank;
2577 #ifdef PASO_MPI
2578 error = MPI_Comm_rank( get_MPIComm(), &rank );
2579 #else
2580 rank = 0;
2581 #endif
2582 return rank;
2583 }
2584
2585 MPI_Comm
2586 Data::get_MPIComm() const
2587 {
2588 #ifdef PASO_MPI
2589 return MPI_COMM_WORLD;
2590 #else
2591 return -1;
2592 #endif
2593 }
2594

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26