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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 150 - (show annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years ago) by jgs
Original Path: trunk/esys2/escript/src/Data/Data.cpp
File size: 50424 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26