/[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 123 - (show annotations)
Fri Jul 8 04:08:13 2005 UTC (13 years, 11 months ago) by jgs
File size: 47581 byte(s)
Merge of development branch back to main trunk on 2005-07-08

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26