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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 122 - (show annotations)
Thu Jun 9 05:38:05 2005 UTC (14 years ago) by jgs
File size: 43332 byte(s)
Merge of development branch back to main trunk on 2005-06-09

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26