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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 751 - (show annotations)
Mon Jun 26 01:46:34 2006 UTC (12 years, 9 months ago) by bcumming
File size: 54964 byte(s)
Changes relating to the MPI version of escript
The standard OpenMP version of escript is unchanged

- updated data types (Finley_Mesh, Finley_NodeFile, etc) to store meshes
  over multiple MPI processes.
- added CommBuffer code in Paso for communication of Data associated
  with distributed meshes
- updates in Finley and Escript to support distributed data and operations
  on distributed data (such as interpolation).
- construction of RHS in MPI, so that simple explicit schemes (such as
  /docs/examples/wave.py without IO and the Locator) can run in MPI.
- updated mesh generation for first order line, rectangle and brick
  meshes and second order line meshes in MPI.        
- small changes to trunk/SConstruct and trunk/scons/ess_options.py to
  build the MPI version, these changes are turned off by default.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26