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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 148 - (show annotations)
Tue Aug 23 01:24:31 2005 UTC (14 years ago) by jgs
Original Path: trunk/esys2/escript/src/Data/Data.cpp
File size: 52180 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-08-23

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26