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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 121 - (show annotations)
Fri May 6 04:26:16 2005 UTC (14 years, 8 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/Data.cpp
File size: 43372 byte(s)
Merge of development branch back to main trunk on 2005-05-06

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26