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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1455 - (hide annotations)
Thu Feb 28 17:19:44 2008 UTC (11 years, 6 months ago) by phornby
File size: 75892 byte(s)
Merge of branches/windows_from_1431_trunk.

Revamp of the exception system.
Fix unused vars and signed/unsigned comparisons.
defined a macro THROW(ARG) in the system_dep.h's to
deal with the expectations of declarations on different architectures.

Details in the logs of branches/windows_from_1431_trunk.

pre-merge snapshot of the trunk in tags/trunk_at_1452


1 jgs 480
2 ksteube 1312 /* $Id$ */
3    
4     /*******************************************************
5     *
6     * Copyright 2003-2007 by ACceSS MNRF
7     * Copyright 2007 by University of Queensland
8     *
9     * http://esscc.uq.edu.au
10     * Primary Business: Queensland, Australia
11     * Licensed under the Open Software License version 3.0
12     * http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16 jgs 474 #include "Data.h"
17 jgs 94
18 jgs 480 #include "DataExpanded.h"
19     #include "DataConstant.h"
20     #include "DataTagged.h"
21     #include "DataEmpty.h"
22     #include "DataArrayView.h"
23     #include "FunctionSpaceFactory.h"
24     #include "AbstractContinuousDomain.h"
25     #include "UnaryFuncs.h"
26 ksteube 1312 extern "C" {
27     #include "escript/blocktimer.h"
28     }
29 jgs 480
30 jgs 119 #include <fstream>
31 jgs 94 #include <algorithm>
32     #include <vector>
33     #include <functional>
34    
35 jgs 153 #include <boost/python/dict.hpp>
36 jgs 94 #include <boost/python/extract.hpp>
37     #include <boost/python/long.hpp>
38    
39     using namespace std;
40     using namespace boost::python;
41     using namespace boost;
42     using namespace escript;
43    
44     Data::Data()
45     {
46     //
47     // Default data is type DataEmpty
48     DataAbstract* temp=new DataEmpty();
49 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
50     m_data=temp_data;
51 gross 783 m_protected=false;
52 jgs 94 }
53    
54     Data::Data(double value,
55     const tuple& shape,
56     const FunctionSpace& what,
57     bool expanded)
58     {
59     DataArrayView::ShapeType dataPointShape;
60     for (int i = 0; i < shape.attr("__len__")(); ++i) {
61     dataPointShape.push_back(extract<const int>(shape[i]));
62     }
63 matt 1319
64     int len = DataArrayView::noValues(dataPointShape);
65     DataVector temp_data(len,value,len);
66     DataArrayView temp_dataView(temp_data, dataPointShape);
67    
68     initialise(temp_dataView, what, expanded);
69    
70 gross 783 m_protected=false;
71 jgs 94 }
72    
73     Data::Data(double value,
74     const DataArrayView::ShapeType& dataPointShape,
75     const FunctionSpace& what,
76     bool expanded)
77     {
78 matt 1319 int len = DataArrayView::noValues(dataPointShape);
79    
80     DataVector temp_data(len,value,len);
81     DataArrayView temp_dataView(temp_data, dataPointShape);
82    
83     initialise(temp_dataView, what, expanded);
84    
85 gross 783 m_protected=false;
86 jgs 94 }
87    
88 jgs 102 Data::Data(const Data& inData)
89 jgs 94 {
90 jgs 102 m_data=inData.m_data;
91 gross 783 m_protected=inData.isProtected();
92 jgs 94 }
93    
94     Data::Data(const Data& inData,
95     const DataArrayView::RegionType& region)
96     {
97     //
98 jgs 102 // Create Data which is a slice of another Data
99     DataAbstract* tmp = inData.m_data->getSlice(region);
100     shared_ptr<DataAbstract> temp_data(tmp);
101     m_data=temp_data;
102 gross 783 m_protected=false;
103 jgs 94 }
104    
105     Data::Data(const Data& inData,
106     const FunctionSpace& functionspace)
107     {
108     if (inData.getFunctionSpace()==functionspace) {
109     m_data=inData.m_data;
110     } else {
111     Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
112 jgs 123 // Note: Must use a reference or pointer to a derived object
113 jgs 94 // in order to get polymorphic behaviour. Shouldn't really
114     // be able to create an instance of AbstractDomain but that was done
115 jgs 123 // as a boost:python work around which may no longer be required.
116 jgs 94 const AbstractDomain& inDataDomain=inData.getDomain();
117     if (inDataDomain==functionspace.getDomain()) {
118     inDataDomain.interpolateOnDomain(tmp,inData);
119     } else {
120     inDataDomain.interpolateACross(tmp,inData);
121     }
122     m_data=tmp.m_data;
123     }
124 gross 783 m_protected=false;
125 jgs 94 }
126    
127     Data::Data(const DataTagged::TagListType& tagKeys,
128     const DataTagged::ValueListType & values,
129     const DataArrayView& defaultValue,
130     const FunctionSpace& what,
131     bool expanded)
132     {
133     DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
134 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
135     m_data=temp_data;
136 gross 783 m_protected=false;
137 jgs 94 if (expanded) {
138     expand();
139     }
140     }
141    
142     Data::Data(const numeric::array& value,
143     const FunctionSpace& what,
144     bool expanded)
145     {
146     initialise(value,what,expanded);
147 gross 783 m_protected=false;
148 jgs 94 }
149    
150     Data::Data(const DataArrayView& value,
151     const FunctionSpace& what,
152     bool expanded)
153     {
154     initialise(value,what,expanded);
155 gross 783 m_protected=false;
156 jgs 94 }
157    
158     Data::Data(const object& value,
159     const FunctionSpace& what,
160     bool expanded)
161     {
162     numeric::array asNumArray(value);
163     initialise(asNumArray,what,expanded);
164 gross 783 m_protected=false;
165 jgs 94 }
166    
167 matt 1319
168 jgs 94 Data::Data(const object& value,
169     const Data& other)
170     {
171 matt 1319
172     numeric::array asNumArray(value);
173    
174    
175     // extract the shape of the numarray
176     DataArrayView::ShapeType tempShape;
177     for (int i=0; i < asNumArray.getrank(); i++) {
178     tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
179     }
180     // get the space for the data vector
181     int len = DataArrayView::noValues(tempShape);
182     DataVector temp_data(len, 0.0, len);
183     DataArrayView temp_dataView(temp_data, tempShape);
184     temp_dataView.copy(asNumArray);
185    
186 jgs 94 //
187     // Create DataConstant using the given value and all other parameters
188     // copied from other. If value is a rank 0 object this Data
189     // will assume the point data shape of other.
190 matt 1319
191     if (temp_dataView.getRank()==0) {
192     int len = DataArrayView::noValues(other.getPointDataView().getShape());
193    
194     DataVector temp2_data(len, temp_dataView(), len);
195 matt 1328 DataArrayView temp2_dataView(temp2_data, other.getPointDataView().getShape());
196 matt 1319 initialise(temp2_dataView, other.getFunctionSpace(), false);
197    
198 jgs 94 } else {
199     //
200     // Create a DataConstant with the same sample shape as other
201 matt 1319 initialise(temp_dataView, other.getFunctionSpace(), false);
202 jgs 94 }
203 gross 783 m_protected=false;
204 jgs 94 }
205    
206 jgs 151 Data::~Data()
207     {
208    
209     }
210    
211 jgs 94 escriptDataC
212     Data::getDataC()
213     {
214     escriptDataC temp;
215     temp.m_dataPtr=(void*)this;
216     return temp;
217     }
218    
219     escriptDataC
220     Data::getDataC() const
221     {
222     escriptDataC temp;
223     temp.m_dataPtr=(void*)this;
224     return temp;
225     }
226    
227 jgs 121 const boost::python::tuple
228 jgs 94 Data::getShapeTuple() const
229     {
230     const DataArrayView::ShapeType& shape=getDataPointShape();
231     switch(getDataPointRank()) {
232     case 0:
233     return make_tuple();
234     case 1:
235     return make_tuple(long_(shape[0]));
236     case 2:
237     return make_tuple(long_(shape[0]),long_(shape[1]));
238     case 3:
239     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]));
240     case 4:
241     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]),long_(shape[3]));
242     default:
243     throw DataException("Error - illegal Data rank.");
244     }
245     }
246     void
247     Data::copy(const Data& other)
248     {
249     //
250     // Perform a deep copy
251     {
252     DataExpanded* temp=dynamic_cast<DataExpanded*>(other.m_data.get());
253     if (temp!=0) {
254     //
255     // Construct a DataExpanded copy
256     DataAbstract* newData=new DataExpanded(*temp);
257 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
258     m_data=temp_data;
259 jgs 94 return;
260     }
261     }
262     {
263     DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());
264     if (temp!=0) {
265     //
266 jgs 102 // Construct a DataTagged copy
267 jgs 94 DataAbstract* newData=new DataTagged(*temp);
268 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
269     m_data=temp_data;
270 jgs 94 return;
271     }
272     }
273     {
274     DataConstant* temp=dynamic_cast<DataConstant*>(other.m_data.get());
275     if (temp!=0) {
276     //
277     // Construct a DataConstant copy
278     DataAbstract* newData=new DataConstant(*temp);
279 jgs 102 shared_ptr<DataAbstract> temp_data(newData);
280     m_data=temp_data;
281 jgs 94 return;
282     }
283     }
284 jgs 102 {
285     DataEmpty* temp=dynamic_cast<DataEmpty*>(other.m_data.get());
286     if (temp!=0) {
287     //
288     // Construct a DataEmpty copy
289     DataAbstract* newData=new DataEmpty();
290     shared_ptr<DataAbstract> temp_data(newData);
291     m_data=temp_data;
292     return;
293     }
294     }
295 jgs 94 throw DataException("Error - Copy not implemented for this Data type.");
296     }
297    
298 gross 1118
299 jgs 94 void
300 gross 1118 Data::setToZero()
301     {
302     {
303     DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
304     if (temp!=0) {
305     temp->setToZero();
306     return;
307     }
308     }
309     {
310     DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
311     if (temp!=0) {
312     temp->setToZero();
313     return;
314     }
315     }
316     {
317     DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
318     if (temp!=0) {
319     temp->setToZero();
320     return;
321     }
322     }
323     throw DataException("Error - Data can not be set to zero.");
324     }
325    
326     void
327 jgs 94 Data::copyWithMask(const Data& other,
328     const Data& mask)
329     {
330     Data mask1;
331     Data mask2;
332    
333     mask1 = mask.wherePositive();
334     mask2.copy(mask1);
335    
336     mask1 *= other;
337     mask2 *= *this;
338     mask2 = *this - mask2;
339    
340     *this = mask1 + mask2;
341     }
342    
343     bool
344     Data::isExpanded() const
345     {
346     DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
347     return (temp!=0);
348     }
349    
350     bool
351     Data::isTagged() const
352     {
353     DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
354     return (temp!=0);
355     }
356    
357     bool
358     Data::isEmpty() const
359     {
360     DataEmpty* temp=dynamic_cast<DataEmpty*>(m_data.get());
361     return (temp!=0);
362     }
363    
364     bool
365     Data::isConstant() const
366     {
367     DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
368     return (temp!=0);
369     }
370    
371     void
372 ksteube 1312 Data::setProtection()
373     {
374 gross 783 m_protected=true;
375     }
376    
377     bool
378 ksteube 1312 Data::isProtected() const
379     {
380 gross 783 return m_protected;
381     }
382    
383    
384    
385     void
386 jgs 94 Data::expand()
387     {
388     if (isConstant()) {
389     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
390     DataAbstract* temp=new DataExpanded(*tempDataConst);
391 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
392     m_data=temp_data;
393 jgs 94 } else if (isTagged()) {
394     DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
395     DataAbstract* temp=new DataExpanded(*tempDataTag);
396 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
397     m_data=temp_data;
398 jgs 94 } else if (isExpanded()) {
399     //
400     // do nothing
401     } else if (isEmpty()) {
402     throw DataException("Error - Expansion of DataEmpty not possible.");
403     } else {
404     throw DataException("Error - Expansion not implemented for this Data type.");
405     }
406     }
407    
408     void
409     Data::tag()
410     {
411     if (isConstant()) {
412     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
413     DataAbstract* temp=new DataTagged(*tempDataConst);
414 jgs 102 shared_ptr<DataAbstract> temp_data(temp);
415     m_data=temp_data;
416 jgs 94 } else if (isTagged()) {
417     // do nothing
418     } else if (isExpanded()) {
419     throw DataException("Error - Creating tag data from DataExpanded not possible.");
420     } else if (isEmpty()) {
421     throw DataException("Error - Creating tag data from DataEmpty not possible.");
422     } else {
423     throw DataException("Error - Tagging not implemented for this Data type.");
424     }
425     }
426    
427 gross 854 Data
428     Data::oneOver() const
429 jgs 102 {
430 matt 1334 return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
431 jgs 102 }
432    
433 jgs 94 Data
434 gross 698 Data::wherePositive() const
435 jgs 94 {
436 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
437 jgs 94 }
438    
439     Data
440 gross 698 Data::whereNegative() const
441 jgs 102 {
442 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
443 jgs 102 }
444    
445     Data
446 gross 698 Data::whereNonNegative() const
447 jgs 94 {
448 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
449 jgs 94 }
450    
451     Data
452 gross 698 Data::whereNonPositive() const
453 jgs 94 {
454 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
455 jgs 94 }
456    
457     Data
458 jgs 571 Data::whereZero(double tol) const
459 jgs 94 {
460 jgs 571 Data dataAbs=abs();
461 matt 1334 return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
462 jgs 94 }
463    
464     Data
465 jgs 571 Data::whereNonZero(double tol) const
466 jgs 102 {
467 jgs 571 Data dataAbs=abs();
468 matt 1334 return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
469 jgs 102 }
470    
471     Data
472 jgs 94 Data::interpolate(const FunctionSpace& functionspace) const
473     {
474     return Data(*this,functionspace);
475     }
476    
477     bool
478     Data::probeInterpolation(const FunctionSpace& functionspace) const
479     {
480     if (getFunctionSpace()==functionspace) {
481     return true;
482     } else {
483     const AbstractDomain& domain=getDomain();
484     if (domain==functionspace.getDomain()) {
485     return domain.probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());
486     } else {
487     return domain.probeInterpolationACross(getFunctionSpace().getTypeCode(),functionspace.getDomain(),functionspace.getTypeCode());
488     }
489     }
490     }
491    
492     Data
493     Data::gradOn(const FunctionSpace& functionspace) const
494     {
495 matt 1353 double blocktimer_start = blocktimer_time();
496 jgs 94 if (functionspace.getDomain()!=getDomain())
497     throw DataException("Error - gradient cannot be calculated on different domains.");
498     DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
499     grad_shape.push_back(functionspace.getDim());
500     Data out(0.0,grad_shape,functionspace,true);
501     getDomain().setToGradient(out,*this);
502 matt 1353 blocktimer_increment("grad()", blocktimer_start);
503 jgs 94 return out;
504     }
505    
506     Data
507     Data::grad() const
508     {
509     return gradOn(escript::function(getDomain()));
510     }
511    
512     int
513     Data::getDataPointSize() const
514     {
515     return getPointDataView().noValues();
516     }
517    
518     DataArrayView::ValueType::size_type
519     Data::getLength() const
520     {
521     return m_data->getLength();
522     }
523    
524     const DataArrayView::ShapeType&
525     Data::getDataPointShape() const
526     {
527     return getPointDataView().getShape();
528     }
529    
530 gross 921
531    
532 ksteube 1312 const
533 jgs 121 boost::python::numeric::array
534 ksteube 1312 Data:: getValueOfDataPoint(int dataPointNo)
535 jgs 121 {
536 gross 921 size_t length=0;
537     int i, j, k, l;
538 jgs 121 //
539     // determine the rank and shape of each data point
540     int dataPointRank = getDataPointRank();
541     DataArrayView::ShapeType dataPointShape = getDataPointShape();
542    
543     //
544     // create the numeric array to be returned
545     boost::python::numeric::array numArray(0.0);
546    
547     //
548 gross 921 // the shape of the returned numeric array will be the same
549     // as that of the data point
550     int arrayRank = dataPointRank;
551     DataArrayView::ShapeType arrayShape = dataPointShape;
552 jgs 121
553     //
554     // resize the numeric array to the shape just calculated
555 gross 921 if (arrayRank==0) {
556     numArray.resize(1);
557     }
558 jgs 121 if (arrayRank==1) {
559     numArray.resize(arrayShape[0]);
560     }
561     if (arrayRank==2) {
562     numArray.resize(arrayShape[0],arrayShape[1]);
563     }
564     if (arrayRank==3) {
565     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
566     }
567     if (arrayRank==4) {
568     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
569     }
570    
571 gross 921 if (getNumDataPointsPerSample()>0) {
572     int sampleNo = dataPointNo/getNumDataPointsPerSample();
573     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
574     //
575     // Check a valid sample number has been supplied
576 trankine 924 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
577 gross 921 throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
578     }
579 ksteube 1312
580 gross 921 //
581     // Check a valid data point number has been supplied
582 trankine 924 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
583 gross 921 throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
584     }
585     // TODO: global error handling
586     // create a view of the data if it is stored locally
587     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
588 ksteube 1312
589 gross 921 switch( dataPointRank ){
590     case 0 :
591     numArray[0] = dataPointView();
592     break;
593 ksteube 1312 case 1 :
594 gross 921 for( i=0; i<dataPointShape[0]; i++ )
595     numArray[i]=dataPointView(i);
596     break;
597 ksteube 1312 case 2 :
598 gross 921 for( i=0; i<dataPointShape[0]; i++ )
599     for( j=0; j<dataPointShape[1]; j++)
600     numArray[make_tuple(i,j)]=dataPointView(i,j);
601     break;
602 ksteube 1312 case 3 :
603 gross 921 for( i=0; i<dataPointShape[0]; i++ )
604     for( j=0; j<dataPointShape[1]; j++ )
605     for( k=0; k<dataPointShape[2]; k++)
606     numArray[make_tuple(i,j,k)]=dataPointView(i,j,k);
607     break;
608     case 4 :
609     for( i=0; i<dataPointShape[0]; i++ )
610     for( j=0; j<dataPointShape[1]; j++ )
611     for( k=0; k<dataPointShape[2]; k++ )
612     for( l=0; l<dataPointShape[3]; l++)
613     numArray[make_tuple(i,j,k,l)]=dataPointView(i,j,k,l);
614     break;
615     }
616 jgs 117 }
617     //
618 gross 921 // return the array
619 jgs 117 return numArray;
620 gross 921
621 jgs 117 }
622 gross 921 void
623 ksteube 1312 Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
624 jgs 121 {
625 gross 1034 // this will throw if the value cannot be represented
626     boost::python::numeric::array num_array(py_object);
627     setValueOfDataPointToArray(dataPointNo,num_array);
628    
629    
630     }
631    
632     void
633     Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)
634     {
635 gross 921 if (isProtected()) {
636     throw DataException("Error - attempt to update protected Data object.");
637     }
638     //
639     // check rank
640 ksteube 1312 if (num_array.getrank()<getDataPointRank())
641 gross 921 throw DataException("Rank of numarray does not match Data object rank");
642 bcumming 790
643 jgs 121 //
644 gross 921 // check shape of num_array
645     for (int i=0; i<getDataPointRank(); i++) {
646     if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
647     throw DataException("Shape of numarray does not match Data object rank");
648 jgs 121 }
649     //
650 gross 921 // make sure data is expanded:
651     if (!isExpanded()) {
652     expand();
653 jgs 121 }
654 gross 921 if (getNumDataPointsPerSample()>0) {
655     int sampleNo = dataPointNo/getNumDataPointsPerSample();
656     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
657     m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
658     } else {
659     m_data->copyToDataPoint(-1, 0,num_array);
660     }
661     }
662 jgs 121
663 gross 922 void
664     Data::setValueOfDataPoint(int dataPointNo, const double value)
665     {
666     if (isProtected()) {
667     throw DataException("Error - attempt to update protected Data object.");
668     }
669     //
670     // make sure data is expanded:
671     if (!isExpanded()) {
672     expand();
673     }
674     if (getNumDataPointsPerSample()>0) {
675     int sampleNo = dataPointNo/getNumDataPointsPerSample();
676     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
677     m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
678     } else {
679     m_data->copyToDataPoint(-1, 0,value);
680     }
681     }
682    
683 ksteube 1312 const
684 gross 921 boost::python::numeric::array
685 ksteube 1312 Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
686 gross 921 {
687     size_t length=0;
688     int i, j, k, l, pos;
689 jgs 121 //
690     // determine the rank and shape of each data point
691     int dataPointRank = getDataPointRank();
692     DataArrayView::ShapeType dataPointShape = getDataPointShape();
693    
694     //
695     // create the numeric array to be returned
696     boost::python::numeric::array numArray(0.0);
697    
698     //
699     // the shape of the returned numeric array will be the same
700     // as that of the data point
701     int arrayRank = dataPointRank;
702     DataArrayView::ShapeType arrayShape = dataPointShape;
703    
704     //
705     // resize the numeric array to the shape just calculated
706     if (arrayRank==0) {
707     numArray.resize(1);
708     }
709     if (arrayRank==1) {
710     numArray.resize(arrayShape[0]);
711     }
712     if (arrayRank==2) {
713     numArray.resize(arrayShape[0],arrayShape[1]);
714     }
715     if (arrayRank==3) {
716     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
717     }
718     if (arrayRank==4) {
719     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
720     }
721    
722 gross 921 // added for the MPI communication
723     length=1;
724     for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
725     double *tmpData = new double[length];
726 bcumming 790
727 jgs 121 //
728     // load the values for the data point into the numeric array.
729 bcumming 790
730     // updated for the MPI case
731     if( get_MPIRank()==procNo ){
732 gross 921 if (getNumDataPointsPerSample()>0) {
733     int sampleNo = dataPointNo/getNumDataPointsPerSample();
734     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
735     //
736     // Check a valid sample number has been supplied
737 trankine 924 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
738 gross 921 throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
739     }
740 ksteube 1312
741 gross 921 //
742     // Check a valid data point number has been supplied
743 trankine 924 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
744 gross 921 throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
745     }
746     // TODO: global error handling
747 bcumming 790 // create a view of the data if it is stored locally
748 gross 921 DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
749 ksteube 1312
750 bcumming 790 // pack the data from the view into tmpData for MPI communication
751     pos=0;
752     switch( dataPointRank ){
753     case 0 :
754     tmpData[0] = dataPointView();
755     break;
756 ksteube 1312 case 1 :
757 bcumming 790 for( i=0; i<dataPointShape[0]; i++ )
758     tmpData[i]=dataPointView(i);
759     break;
760 ksteube 1312 case 2 :
761 bcumming 790 for( i=0; i<dataPointShape[0]; i++ )
762     for( j=0; j<dataPointShape[1]; j++, pos++ )
763     tmpData[pos]=dataPointView(i,j);
764     break;
765 ksteube 1312 case 3 :
766 bcumming 790 for( i=0; i<dataPointShape[0]; i++ )
767     for( j=0; j<dataPointShape[1]; j++ )
768     for( k=0; k<dataPointShape[2]; k++, pos++ )
769     tmpData[pos]=dataPointView(i,j,k);
770     break;
771     case 4 :
772     for( i=0; i<dataPointShape[0]; i++ )
773     for( j=0; j<dataPointShape[1]; j++ )
774     for( k=0; k<dataPointShape[2]; k++ )
775     for( l=0; l<dataPointShape[3]; l++, pos++ )
776     tmpData[pos]=dataPointView(i,j,k,l);
777     break;
778     }
779 gross 921 }
780 bcumming 790 }
781 ksteube 1312 #ifdef PASO_MPI
782 gross 921 // broadcast the data to all other processes
783     MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
784     #endif
785 bcumming 790
786     // unpack the data
787     switch( dataPointRank ){
788     case 0 :
789 gross 921 numArray[0]=tmpData[0];
790 bcumming 790 break;
791 ksteube 1312 case 1 :
792 bcumming 790 for( i=0; i<dataPointShape[0]; i++ )
793     numArray[i]=tmpData[i];
794     break;
795 ksteube 1312 case 2 :
796 bcumming 790 for( i=0; i<dataPointShape[0]; i++ )
797     for( j=0; j<dataPointShape[1]; j++ )
798 gross 921 numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
799 bcumming 790 break;
800 ksteube 1312 case 3 :
801 bcumming 790 for( i=0; i<dataPointShape[0]; i++ )
802     for( j=0; j<dataPointShape[1]; j++ )
803     for( k=0; k<dataPointShape[2]; k++ )
804 gross 921 numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
805 bcumming 790 break;
806     case 4 :
807     for( i=0; i<dataPointShape[0]; i++ )
808     for( j=0; j<dataPointShape[1]; j++ )
809     for( k=0; k<dataPointShape[2]; k++ )
810     for( l=0; l<dataPointShape[3]; l++ )
811 gross 921 numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
812 bcumming 790 break;
813     }
814    
815 ksteube 1312 delete [] tmpData;
816 jgs 121 //
817     // return the loaded array
818     return numArray;
819     }
820    
821 gross 921
822    
823 jgs 121 boost::python::numeric::array
824 jgs 94 Data::integrate() const
825     {
826     int index;
827     int rank = getDataPointRank();
828     DataArrayView::ShapeType shape = getDataPointShape();
829 ksteube 1312 int dataPointSize = getDataPointSize();
830 jgs 94
831     //
832     // calculate the integral values
833 ksteube 1312 vector<double> integrals(dataPointSize);
834     vector<double> integrals_local(dataPointSize);
835     #ifdef PASO_MPI
836     AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals_local,*this);
837     // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
838     double *tmp = new double[dataPointSize];
839     double *tmp_local = new double[dataPointSize];
840     for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
841     MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
842     for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
843     delete[] tmp;
844     delete[] tmp_local;
845     #else
846 jgs 94 AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);
847 ksteube 1312 #endif
848 jgs 94
849     //
850     // create the numeric array to be returned
851     // and load the array with the integral values
852     boost::python::numeric::array bp_array(1.0);
853     if (rank==0) {
854 jgs 108 bp_array.resize(1);
855 jgs 94 index = 0;
856     bp_array[0] = integrals[index];
857     }
858     if (rank==1) {
859     bp_array.resize(shape[0]);
860     for (int i=0; i<shape[0]; i++) {
861     index = i;
862     bp_array[i] = integrals[index];
863     }
864     }
865     if (rank==2) {
866 gross 436 bp_array.resize(shape[0],shape[1]);
867     for (int i=0; i<shape[0]; i++) {
868     for (int j=0; j<shape[1]; j++) {
869     index = i + shape[0] * j;
870     bp_array[make_tuple(i,j)] = integrals[index];
871     }
872     }
873 jgs 94 }
874     if (rank==3) {
875     bp_array.resize(shape[0],shape[1],shape[2]);
876     for (int i=0; i<shape[0]; i++) {
877     for (int j=0; j<shape[1]; j++) {
878     for (int k=0; k<shape[2]; k++) {
879     index = i + shape[0] * ( j + shape[1] * k );
880 gross 436 bp_array[make_tuple(i,j,k)] = integrals[index];
881 jgs 94 }
882     }
883     }
884     }
885     if (rank==4) {
886     bp_array.resize(shape[0],shape[1],shape[2],shape[3]);
887     for (int i=0; i<shape[0]; i++) {
888     for (int j=0; j<shape[1]; j++) {
889     for (int k=0; k<shape[2]; k++) {
890     for (int l=0; l<shape[3]; l++) {
891     index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
892 gross 436 bp_array[make_tuple(i,j,k,l)] = integrals[index];
893 jgs 94 }
894     }
895     }
896     }
897     }
898    
899     //
900     // return the loaded array
901     return bp_array;
902     }
903    
904     Data
905     Data::sin() const
906     {
907 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
908 jgs 94 }
909    
910     Data
911     Data::cos() const
912     {
913 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
914 jgs 94 }
915    
916     Data
917     Data::tan() const
918     {
919 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
920 jgs 94 }
921    
922     Data
923 jgs 150 Data::asin() const
924     {
925 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
926 jgs 150 }
927    
928     Data
929     Data::acos() const
930     {
931 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
932 jgs 150 }
933    
934 phornby 1026
935 jgs 150 Data
936     Data::atan() const
937     {
938 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
939 jgs 150 }
940    
941     Data
942     Data::sinh() const
943     {
944 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
945 matt 1334
946 jgs 150 }
947    
948     Data
949     Data::cosh() const
950     {
951 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
952 jgs 150 }
953    
954     Data
955     Data::tanh() const
956     {
957 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
958 jgs 150 }
959    
960 phornby 1026
961 jgs 150 Data
962 phornby 1026 Data::erf() const
963     {
964 gross 1028 #ifdef _WIN32
965     throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
966     #else
967 matt 1334 return C_TensorUnaryOperation(*this, ::erf);
968 phornby 1026 #endif
969     }
970    
971     Data
972 jgs 150 Data::asinh() const
973     {
974 phornby 1032 #ifdef _WIN32
975 matt 1334 return C_TensorUnaryOperation(*this, escript::asinh_substitute);
976 phornby 1032 #else
977 matt 1334 return C_TensorUnaryOperation(*this, ::asinh);
978 phornby 1032 #endif
979 jgs 150 }
980    
981     Data
982     Data::acosh() const
983     {
984 phornby 1032 #ifdef _WIN32
985 matt 1334 return C_TensorUnaryOperation(*this, escript::acosh_substitute);
986 phornby 1032 #else
987 matt 1334 return C_TensorUnaryOperation(*this, ::acosh);
988 phornby 1032 #endif
989 jgs 150 }
990    
991     Data
992     Data::atanh() const
993     {
994 phornby 1032 #ifdef _WIN32
995 matt 1334 return C_TensorUnaryOperation(*this, escript::atanh_substitute);
996 phornby 1032 #else
997 matt 1334 return C_TensorUnaryOperation(*this, ::atanh);
998 phornby 1032 #endif
999 jgs 150 }
1000    
1001     Data
1002 gross 286 Data::log10() const
1003 jgs 94 {
1004 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1005 jgs 94 }
1006    
1007     Data
1008 gross 286 Data::log() const
1009 jgs 94 {
1010 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1011 jgs 94 }
1012    
1013 jgs 106 Data
1014     Data::sign() const
1015 jgs 94 {
1016 matt 1334 return C_TensorUnaryOperation(*this, escript::fsign);
1017 jgs 94 }
1018    
1019 jgs 106 Data
1020     Data::abs() const
1021 jgs 94 {
1022 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1023 jgs 94 }
1024    
1025 jgs 106 Data
1026     Data::neg() const
1027 jgs 94 {
1028 matt 1334 return C_TensorUnaryOperation(*this, negate<double>());
1029 jgs 94 }
1030    
1031 jgs 102 Data
1032 jgs 106 Data::pos() const
1033 jgs 94 {
1034 jgs 148 Data result;
1035     // perform a deep copy
1036     result.copy(*this);
1037     return result;
1038 jgs 102 }
1039    
1040     Data
1041 jgs 106 Data::exp() const
1042 jgs 102 {
1043 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1044 jgs 102 }
1045    
1046     Data
1047 jgs 106 Data::sqrt() const
1048 jgs 102 {
1049 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1050 jgs 102 }
1051    
1052 jgs 106 double
1053     Data::Lsup() const
1054 jgs 102 {
1055 phornby 1455 double localValue;
1056 jgs 106 //
1057     // set the initial absolute maximum value to zero
1058 bcumming 751
1059 jgs 147 AbsMax abs_max_func;
1060 bcumming 751 localValue = algorithm(abs_max_func,0);
1061     #ifdef PASO_MPI
1062 phornby 1455 double globalValue;
1063 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1064     return globalValue;
1065     #else
1066     return localValue;
1067     #endif
1068 jgs 117 }
1069    
1070     double
1071 jgs 106 Data::sup() const
1072 jgs 102 {
1073 phornby 1455 double localValue;
1074 jgs 106 //
1075     // set the initial maximum value to min possible double
1076 jgs 147 FMax fmax_func;
1077 bcumming 751 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1078     #ifdef PASO_MPI
1079 phornby 1455 double globalValue;
1080 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1081     return globalValue;
1082     #else
1083     return localValue;
1084     #endif
1085 jgs 102 }
1086    
1087 jgs 106 double
1088     Data::inf() const
1089 jgs 102 {
1090 phornby 1455 double localValue;
1091 jgs 106 //
1092     // set the initial minimum value to max possible double
1093 jgs 147 FMin fmin_func;
1094 bcumming 751 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1095     #ifdef PASO_MPI
1096 phornby 1455 double globalValue;
1097 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1098     return globalValue;
1099     #else
1100     return localValue;
1101     #endif
1102 jgs 102 }
1103    
1104 bcumming 751 /* TODO */
1105     /* global reduction */
1106 jgs 102 Data
1107 jgs 106 Data::maxval() const
1108 jgs 102 {
1109 jgs 113 //
1110     // set the initial maximum value to min possible double
1111 jgs 147 FMax fmax_func;
1112     return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1113 jgs 102 }
1114    
1115     Data
1116 jgs 106 Data::minval() const
1117 jgs 102 {
1118 jgs 113 //
1119     // set the initial minimum value to max possible double
1120 jgs 147 FMin fmin_func;
1121     return dp_algorithm(fmin_func,numeric_limits<double>::max());
1122 jgs 102 }
1123    
1124 jgs 123 Data
1125 gross 804 Data::swapaxes(const int axis0, const int axis1) const
1126 jgs 123 {
1127 gross 804 int axis0_tmp,axis1_tmp;
1128 gross 800 DataArrayView::ShapeType s=getDataPointShape();
1129     DataArrayView::ShapeType ev_shape;
1130     // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1131     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1132     int rank=getDataPointRank();
1133 gross 804 if (rank<2) {
1134     throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1135 gross 800 }
1136 gross 804 if (axis0<0 || axis0>rank-1) {
1137     throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1138     }
1139     if (axis1<0 || axis1>rank-1) {
1140     throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1141     }
1142     if (axis0 == axis1) {
1143     throw DataException("Error - Data::swapaxes: axis indices must be different.");
1144     }
1145     if (axis0 > axis1) {
1146     axis0_tmp=axis1;
1147     axis1_tmp=axis0;
1148     } else {
1149     axis0_tmp=axis0;
1150     axis1_tmp=axis1;
1151     }
1152 gross 800 for (int i=0; i<rank; i++) {
1153 gross 804 if (i == axis0_tmp) {
1154 ksteube 1312 ev_shape.push_back(s[axis1_tmp]);
1155 gross 804 } else if (i == axis1_tmp) {
1156 ksteube 1312 ev_shape.push_back(s[axis0_tmp]);
1157 gross 800 } else {
1158 ksteube 1312 ev_shape.push_back(s[i]);
1159 gross 800 }
1160     }
1161     Data ev(0.,ev_shape,getFunctionSpace());
1162     ev.typeMatchRight(*this);
1163 gross 804 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1164 gross 800 return ev;
1165    
1166 jgs 123 }
1167    
1168     Data
1169 ksteube 775 Data::symmetric() const
1170 jgs 123 {
1171 ksteube 775 // check input
1172     DataArrayView::ShapeType s=getDataPointShape();
1173     if (getDataPointRank()==2) {
1174 ksteube 1312 if(s[0] != s[1])
1175 ksteube 775 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1176     }
1177     else if (getDataPointRank()==4) {
1178     if(!(s[0] == s[2] && s[1] == s[3]))
1179     throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1180     }
1181     else {
1182     throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1183     }
1184     Data ev(0.,getDataPointShape(),getFunctionSpace());
1185     ev.typeMatchRight(*this);
1186     m_data->symmetric(ev.m_data.get());
1187     return ev;
1188     }
1189    
1190     Data
1191     Data::nonsymmetric() const
1192     {
1193     // check input
1194     DataArrayView::ShapeType s=getDataPointShape();
1195     if (getDataPointRank()==2) {
1196 ksteube 1312 if(s[0] != s[1])
1197 ksteube 775 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1198     DataArrayView::ShapeType ev_shape;
1199     ev_shape.push_back(s[0]);
1200     ev_shape.push_back(s[1]);
1201     Data ev(0.,ev_shape,getFunctionSpace());
1202     ev.typeMatchRight(*this);
1203     m_data->nonsymmetric(ev.m_data.get());
1204     return ev;
1205     }
1206     else if (getDataPointRank()==4) {
1207     if(!(s[0] == s[2] && s[1] == s[3]))
1208     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1209     DataArrayView::ShapeType ev_shape;
1210     ev_shape.push_back(s[0]);
1211     ev_shape.push_back(s[1]);
1212     ev_shape.push_back(s[2]);
1213     ev_shape.push_back(s[3]);
1214     Data ev(0.,ev_shape,getFunctionSpace());
1215     ev.typeMatchRight(*this);
1216     m_data->nonsymmetric(ev.m_data.get());
1217     return ev;
1218     }
1219     else {
1220     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1221     }
1222     }
1223    
1224     Data
1225 gross 800 Data::trace(int axis_offset) const
1226 ksteube 775 {
1227     DataArrayView::ShapeType s=getDataPointShape();
1228     if (getDataPointRank()==2) {
1229     DataArrayView::ShapeType ev_shape;
1230     Data ev(0.,ev_shape,getFunctionSpace());
1231     ev.typeMatchRight(*this);
1232 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1233 ksteube 775 return ev;
1234     }
1235     if (getDataPointRank()==3) {
1236     DataArrayView::ShapeType ev_shape;
1237     if (axis_offset==0) {
1238     int s2=s[2];
1239     ev_shape.push_back(s2);
1240     }
1241     else if (axis_offset==1) {
1242     int s0=s[0];
1243     ev_shape.push_back(s0);
1244     }
1245     Data ev(0.,ev_shape,getFunctionSpace());
1246     ev.typeMatchRight(*this);
1247 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1248 ksteube 775 return ev;
1249     }
1250     if (getDataPointRank()==4) {
1251     DataArrayView::ShapeType ev_shape;
1252     if (axis_offset==0) {
1253     ev_shape.push_back(s[2]);
1254     ev_shape.push_back(s[3]);
1255     }
1256     else if (axis_offset==1) {
1257     ev_shape.push_back(s[0]);
1258     ev_shape.push_back(s[3]);
1259     }
1260     else if (axis_offset==2) {
1261     ev_shape.push_back(s[0]);
1262     ev_shape.push_back(s[1]);
1263     }
1264     Data ev(0.,ev_shape,getFunctionSpace());
1265     ev.typeMatchRight(*this);
1266 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1267 ksteube 775 return ev;
1268     }
1269     else {
1270 gross 800 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1271 ksteube 775 }
1272     }
1273    
1274     Data
1275     Data::transpose(int axis_offset) const
1276     {
1277     DataArrayView::ShapeType s=getDataPointShape();
1278     DataArrayView::ShapeType ev_shape;
1279     // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1280     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1281     int rank=getDataPointRank();
1282     if (axis_offset<0 || axis_offset>rank) {
1283     throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1284     }
1285     for (int i=0; i<rank; i++) {
1286     int index = (axis_offset+i)%rank;
1287     ev_shape.push_back(s[index]); // Append to new shape
1288     }
1289     Data ev(0.,ev_shape,getFunctionSpace());
1290     ev.typeMatchRight(*this);
1291     m_data->transpose(ev.m_data.get(), axis_offset);
1292     return ev;
1293 jgs 123 }
1294    
1295 gross 576 Data
1296     Data::eigenvalues() const
1297     {
1298     // check input
1299     DataArrayView::ShapeType s=getDataPointShape();
1300 ksteube 1312 if (getDataPointRank()!=2)
1301 gross 576 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1302 ksteube 1312 if(s[0] != s[1])
1303 gross 576 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1304     // create return
1305     DataArrayView::ShapeType ev_shape(1,s[0]);
1306     Data ev(0.,ev_shape,getFunctionSpace());
1307     ev.typeMatchRight(*this);
1308     m_data->eigenvalues(ev.m_data.get());
1309     return ev;
1310     }
1311    
1312 jgs 121 const boost::python::tuple
1313 gross 576 Data::eigenvalues_and_eigenvectors(const double tol) const
1314     {
1315     DataArrayView::ShapeType s=getDataPointShape();
1316 ksteube 1312 if (getDataPointRank()!=2)
1317 gross 576 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1318 ksteube 1312 if(s[0] != s[1])
1319 gross 576 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1320     // create return
1321     DataArrayView::ShapeType ev_shape(1,s[0]);
1322     Data ev(0.,ev_shape,getFunctionSpace());
1323     ev.typeMatchRight(*this);
1324     DataArrayView::ShapeType V_shape(2,s[0]);
1325     Data V(0.,V_shape,getFunctionSpace());
1326     V.typeMatchRight(*this);
1327     m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1328     return make_tuple(boost::python::object(ev),boost::python::object(V));
1329     }
1330    
1331     const boost::python::tuple
1332 gross 921 Data::minGlobalDataPoint() const
1333 jgs 121 {
1334 gross 921 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1335 jgs 148 // abort (for unknown reasons) if there are openmp directives with it in the
1336     // surrounding function
1337    
1338     int DataPointNo;
1339 gross 921 int ProcNo;
1340     calc_minGlobalDataPoint(ProcNo,DataPointNo);
1341     return make_tuple(ProcNo,DataPointNo);
1342 jgs 148 }
1343    
1344     void
1345 gross 921 Data::calc_minGlobalDataPoint(int& ProcNo,
1346     int& DataPointNo) const
1347 jgs 148 {
1348     int i,j;
1349     int lowi=0,lowj=0;
1350     double min=numeric_limits<double>::max();
1351    
1352 jgs 121 Data temp=minval();
1353    
1354     int numSamples=temp.getNumSamples();
1355     int numDPPSample=temp.getNumDataPointsPerSample();
1356    
1357 jgs 148 double next,local_min;
1358     int local_lowi,local_lowj;
1359 jgs 121
1360 jgs 148 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1361     {
1362     local_min=min;
1363     #pragma omp for private(i,j) schedule(static)
1364     for (i=0; i<numSamples; i++) {
1365     for (j=0; j<numDPPSample; j++) {
1366     next=temp.getDataPoint(i,j)();
1367     if (next<local_min) {
1368     local_min=next;
1369     local_lowi=i;
1370     local_lowj=j;
1371     }
1372 jgs 121 }
1373     }
1374 jgs 148 #pragma omp critical
1375     if (local_min<min) {
1376     min=local_min;
1377     lowi=local_lowi;
1378     lowj=local_lowj;
1379     }
1380 jgs 121 }
1381    
1382 bcumming 782 #ifdef PASO_MPI
1383     // determine the processor on which the minimum occurs
1384     next = temp.getDataPoint(lowi,lowj)();
1385     int lowProc = 0;
1386     double *globalMins = new double[get_MPISize()+1];
1387     int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1388 ksteube 1312
1389 bcumming 782 if( get_MPIRank()==0 ){
1390     next = globalMins[lowProc];
1391     for( i=1; i<get_MPISize(); i++ )
1392     if( next>globalMins[i] ){
1393     lowProc = i;
1394     next = globalMins[i];
1395     }
1396     }
1397     MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1398    
1399     delete [] globalMins;
1400     ProcNo = lowProc;
1401 bcumming 790 #else
1402     ProcNo = 0;
1403 bcumming 782 #endif
1404 gross 921 DataPointNo = lowj + lowi * numDPPSample;
1405 jgs 121 }
1406    
1407 jgs 104 void
1408     Data::saveDX(std::string fileName) const
1409     {
1410 jgs 153 boost::python::dict args;
1411     args["data"]=boost::python::object(this);
1412     getDomain().saveDX(fileName,args);
1413 jgs 104 return;
1414     }
1415    
1416 jgs 110 void
1417     Data::saveVTK(std::string fileName) const
1418     {
1419 jgs 153 boost::python::dict args;
1420     args["data"]=boost::python::object(this);
1421     getDomain().saveVTK(fileName,args);
1422 jgs 110 return;
1423     }
1424    
1425 jgs 102 Data&
1426     Data::operator+=(const Data& right)
1427     {
1428 gross 783 if (isProtected()) {
1429     throw DataException("Error - attempt to update protected Data object.");
1430     }
1431 jgs 94 binaryOp(right,plus<double>());
1432     return (*this);
1433     }
1434    
1435 jgs 102 Data&
1436     Data::operator+=(const boost::python::object& right)
1437 jgs 94 {
1438 gross 854 Data tmp(right,getFunctionSpace(),false);
1439     binaryOp(tmp,plus<double>());
1440 jgs 94 return (*this);
1441     }
1442 ksteube 1312 Data&
1443     Data::operator=(const Data& other)
1444     {
1445     copy(other);
1446     return (*this);
1447     }
1448 jgs 94
1449 jgs 102 Data&
1450     Data::operator-=(const Data& right)
1451 jgs 94 {
1452 gross 783 if (isProtected()) {
1453     throw DataException("Error - attempt to update protected Data object.");
1454     }
1455 jgs 94 binaryOp(right,minus<double>());
1456     return (*this);
1457     }
1458    
1459 jgs 102 Data&
1460     Data::operator-=(const boost::python::object& right)
1461 jgs 94 {
1462 gross 854 Data tmp(right,getFunctionSpace(),false);
1463     binaryOp(tmp,minus<double>());
1464 jgs 94 return (*this);
1465     }
1466    
1467 jgs 102 Data&
1468     Data::operator*=(const Data& right)
1469 jgs 94 {
1470 gross 783 if (isProtected()) {
1471     throw DataException("Error - attempt to update protected Data object.");
1472     }
1473 jgs 94 binaryOp(right,multiplies<double>());
1474     return (*this);
1475     }
1476    
1477 jgs 102 Data&
1478     Data::operator*=(const boost::python::object& right)
1479 jgs 94 {
1480 gross 854 Data tmp(right,getFunctionSpace(),false);
1481     binaryOp(tmp,multiplies<double>());
1482 jgs 94 return (*this);
1483     }
1484    
1485 jgs 102 Data&
1486     Data::operator/=(const Data& right)
1487 jgs 94 {
1488 gross 783 if (isProtected()) {
1489     throw DataException("Error - attempt to update protected Data object.");
1490     }
1491 jgs 94 binaryOp(right,divides<double>());
1492     return (*this);
1493     }
1494    
1495 jgs 102 Data&
1496     Data::operator/=(const boost::python::object& right)
1497 jgs 94 {
1498 gross 854 Data tmp(right,getFunctionSpace(),false);
1499     binaryOp(tmp,divides<double>());
1500 jgs 94 return (*this);
1501     }
1502    
1503 jgs 102 Data
1504 gross 699 Data::rpowO(const boost::python::object& left) const
1505     {
1506     Data left_d(left,*this);
1507     return left_d.powD(*this);
1508     }
1509    
1510     Data
1511 jgs 102 Data::powO(const boost::python::object& right) const
1512 jgs 94 {
1513 gross 854 Data tmp(right,getFunctionSpace(),false);
1514     return powD(tmp);
1515 jgs 94 }
1516    
1517 jgs 102 Data
1518     Data::powD(const Data& right) const
1519 jgs 94 {
1520 matt 1350 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
1521 jgs 94 }
1522    
1523     //
1524 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1525 jgs 102 Data
1526     escript::operator+(const Data& left, const Data& right)
1527 jgs 94 {
1528 matt 1327 return C_TensorBinaryOperation(left, right, plus<double>());
1529 jgs 94 }
1530    
1531     //
1532 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1533 jgs 102 Data
1534     escript::operator-(const Data& left, const Data& right)
1535 jgs 94 {
1536 matt 1327 return C_TensorBinaryOperation(left, right, minus<double>());
1537 jgs 94 }
1538    
1539     //
1540 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1541 jgs 102 Data
1542     escript::operator*(const Data& left, const Data& right)
1543 jgs 94 {
1544 matt 1327 return C_TensorBinaryOperation(left, right, multiplies<double>());
1545 jgs 94 }
1546    
1547     //
1548 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1549 jgs 102 Data
1550     escript::operator/(const Data& left, const Data& right)
1551 jgs 94 {
1552 matt 1327 return C_TensorBinaryOperation(left, right, divides<double>());
1553 jgs 94 }
1554    
1555     //
1556 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1557 jgs 102 Data
1558     escript::operator+(const Data& left, const boost::python::object& right)
1559 jgs 94 {
1560 gross 854 return left+Data(right,left.getFunctionSpace(),false);
1561 jgs 94 }
1562    
1563     //
1564 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1565 jgs 102 Data
1566     escript::operator-(const Data& left, const boost::python::object& right)
1567 jgs 94 {
1568 gross 854 return left-Data(right,left.getFunctionSpace(),false);
1569 jgs 94 }
1570    
1571     //
1572 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1573 jgs 102 Data
1574     escript::operator*(const Data& left, const boost::python::object& right)
1575 jgs 94 {
1576 gross 854 return left*Data(right,left.getFunctionSpace(),false);
1577 jgs 94 }
1578    
1579     //
1580 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1581 jgs 102 Data
1582     escript::operator/(const Data& left, const boost::python::object& right)
1583 jgs 94 {
1584 gross 854 return left/Data(right,left.getFunctionSpace(),false);
1585 jgs 94 }
1586    
1587     //
1588 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1589 jgs 102 Data
1590     escript::operator+(const boost::python::object& left, const Data& right)
1591 jgs 94 {
1592 gross 854 return Data(left,right.getFunctionSpace(),false)+right;
1593 jgs 94 }
1594    
1595     //
1596 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1597 jgs 102 Data
1598     escript::operator-(const boost::python::object& left, const Data& right)
1599 jgs 94 {
1600 gross 854 return Data(left,right.getFunctionSpace(),false)-right;
1601 jgs 94 }
1602    
1603     //
1604 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1605 jgs 102 Data
1606     escript::operator*(const boost::python::object& left, const Data& right)
1607 jgs 94 {
1608 gross 854 return Data(left,right.getFunctionSpace(),false)*right;
1609 jgs 94 }
1610    
1611     //
1612 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1613 jgs 102 Data
1614     escript::operator/(const boost::python::object& left, const Data& right)
1615 jgs 94 {
1616 gross 854 return Data(left,right.getFunctionSpace(),false)/right;
1617 jgs 94 }
1618    
1619     //
1620 jgs 102 //bool escript::operator==(const Data& left, const Data& right)
1621     //{
1622     // /*
1623     // NB: this operator does very little at this point, and isn't to
1624     // be relied on. Requires further implementation.
1625     // */
1626     //
1627     // bool ret;
1628     //
1629     // if (left.isEmpty()) {
1630     // if(!right.isEmpty()) {
1631     // ret = false;
1632     // } else {
1633     // ret = true;
1634     // }
1635     // }
1636     //
1637     // if (left.isConstant()) {
1638     // if(!right.isConstant()) {
1639     // ret = false;
1640     // } else {
1641     // ret = true;
1642     // }
1643     // }
1644     //
1645     // if (left.isTagged()) {
1646     // if(!right.isTagged()) {
1647     // ret = false;
1648     // } else {
1649     // ret = true;
1650     // }
1651     // }
1652     //
1653     // if (left.isExpanded()) {
1654     // if(!right.isExpanded()) {
1655     // ret = false;
1656     // } else {
1657     // ret = true;
1658     // }
1659     // }
1660     //
1661     // return ret;
1662     //}
1663    
1664 bcumming 751 /* TODO */
1665     /* global reduction */
1666 jgs 102 Data
1667 ksteube 1312 Data::getItem(const boost::python::object& key) const
1668 jgs 94 {
1669 jgs 102 const DataArrayView& view=getPointDataView();
1670 jgs 94
1671 jgs 102 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1672 jgs 94
1673 jgs 102 if (slice_region.size()!=view.getRank()) {
1674     throw DataException("Error - slice size does not match Data rank.");
1675 jgs 94 }
1676    
1677 jgs 102 return getSlice(slice_region);
1678 jgs 94 }
1679    
1680 bcumming 751 /* TODO */
1681     /* global reduction */
1682 jgs 94 Data
1683 jgs 102 Data::getSlice(const DataArrayView::RegionType& region) const
1684 jgs 94 {
1685 jgs 102 return Data(*this,region);
1686 jgs 94 }
1687    
1688 bcumming 751 /* TODO */
1689     /* global reduction */
1690 jgs 94 void
1691 jgs 102 Data::setItemO(const boost::python::object& key,
1692     const boost::python::object& value)
1693 jgs 94 {
1694 jgs 102 Data tempData(value,getFunctionSpace());
1695     setItemD(key,tempData);
1696     }
1697    
1698     void
1699     Data::setItemD(const boost::python::object& key,
1700     const Data& value)
1701     {
1702 jgs 94 const DataArrayView& view=getPointDataView();
1703 jgs 123
1704 jgs 94 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1705     if (slice_region.size()!=view.getRank()) {
1706     throw DataException("Error - slice size does not match Data rank.");
1707     }
1708 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
1709     setSlice(Data(value,getFunctionSpace()),slice_region);
1710     } else {
1711     setSlice(value,slice_region);
1712     }
1713 jgs 94 }
1714    
1715     void
1716 jgs 102 Data::setSlice(const Data& value,
1717     const DataArrayView::RegionType& region)
1718 jgs 94 {
1719 gross 783 if (isProtected()) {
1720     throw DataException("Error - attempt to update protected Data object.");
1721     }
1722 jgs 102 Data tempValue(value);
1723     typeMatchLeft(tempValue);
1724     typeMatchRight(tempValue);
1725     m_data->setSlice(tempValue.m_data.get(),region);
1726     }
1727    
1728     void
1729     Data::typeMatchLeft(Data& right) const
1730     {
1731     if (isExpanded()){
1732     right.expand();
1733     } else if (isTagged()) {
1734     if (right.isConstant()) {
1735     right.tag();
1736     }
1737     }
1738     }
1739    
1740     void
1741     Data::typeMatchRight(const Data& right)
1742     {
1743 jgs 94 if (isTagged()) {
1744     if (right.isExpanded()) {
1745     expand();
1746     }
1747     } else if (isConstant()) {
1748     if (right.isExpanded()) {
1749     expand();
1750     } else if (right.isTagged()) {
1751     tag();
1752     }
1753     }
1754     }
1755    
1756     void
1757 gross 1044 Data::setTaggedValueByName(std::string name,
1758 ksteube 1312 const boost::python::object& value)
1759 gross 1044 {
1760     if (getFunctionSpace().getDomain().isValidTagName(name)) {
1761     int tagKey=getFunctionSpace().getDomain().getTag(name);
1762     setTaggedValue(tagKey,value);
1763     }
1764     }
1765     void
1766 jgs 94 Data::setTaggedValue(int tagKey,
1767     const boost::python::object& value)
1768     {
1769 gross 783 if (isProtected()) {
1770     throw DataException("Error - attempt to update protected Data object.");
1771     }
1772 jgs 94 //
1773     // Ensure underlying data object is of type DataTagged
1774 gross 1358 if (isConstant()) tag();
1775 jgs 94
1776 matt 1319 numeric::array asNumArray(value);
1777 jgs 94
1778 matt 1319
1779     // extract the shape of the numarray
1780     DataArrayView::ShapeType tempShape;
1781     for (int i=0; i < asNumArray.getrank(); i++) {
1782     tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
1783     }
1784    
1785     // get the space for the data vector
1786     int len = DataArrayView::noValues(tempShape);
1787     DataVector temp_data(len, 0.0, len);
1788     DataArrayView temp_dataView(temp_data, tempShape);
1789     temp_dataView.copy(asNumArray);
1790    
1791 jgs 94 //
1792     // Call DataAbstract::setTaggedValue
1793 matt 1319 m_data->setTaggedValue(tagKey,temp_dataView);
1794 jgs 94 }
1795    
1796 jgs 110 void
1797 jgs 121 Data::setTaggedValueFromCPP(int tagKey,
1798     const DataArrayView& value)
1799     {
1800 gross 783 if (isProtected()) {
1801     throw DataException("Error - attempt to update protected Data object.");
1802     }
1803 jgs 121 //
1804     // Ensure underlying data object is of type DataTagged
1805 gross 1358 if (isConstant()) tag();
1806 jgs 121
1807     //
1808     // Call DataAbstract::setTaggedValue
1809     m_data->setTaggedValue(tagKey,value);
1810     }
1811    
1812 jgs 149 int
1813     Data::getTagNumber(int dpno)
1814     {
1815 gross 1409 return getFunctionSpace().getTagFromDataPointNo(dpno);
1816 jgs 149 }
1817    
1818 jgs 121 void
1819 jgs 119 Data::archiveData(const std::string fileName)
1820     {
1821     cout << "Archiving Data object to: " << fileName << endl;
1822    
1823     //
1824     // Determine type of this Data object
1825     int dataType = -1;
1826    
1827     if (isEmpty()) {
1828     dataType = 0;
1829     cout << "\tdataType: DataEmpty" << endl;
1830     }
1831     if (isConstant()) {
1832     dataType = 1;
1833     cout << "\tdataType: DataConstant" << endl;
1834     }
1835     if (isTagged()) {
1836     dataType = 2;
1837     cout << "\tdataType: DataTagged" << endl;
1838     }
1839     if (isExpanded()) {
1840     dataType = 3;
1841     cout << "\tdataType: DataExpanded" << endl;
1842     }
1843 jgs 123
1844 jgs 119 if (dataType == -1) {
1845     throw DataException("archiveData Error: undefined dataType");
1846     }
1847    
1848     //
1849     // Collect data items common to all Data types
1850     int noSamples = getNumSamples();
1851     int noDPPSample = getNumDataPointsPerSample();
1852     int functionSpaceType = getFunctionSpace().getTypeCode();
1853     int dataPointRank = getDataPointRank();
1854     int dataPointSize = getDataPointSize();
1855     int dataLength = getLength();
1856     DataArrayView::ShapeType dataPointShape = getDataPointShape();
1857 woo409 757 vector<int> referenceNumbers(noSamples);
1858 jgs 119 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1859 gross 964 referenceNumbers[sampleNo] = getFunctionSpace().getReferenceIDOfSample(sampleNo);
1860 jgs 119 }
1861 woo409 757 vector<int> tagNumbers(noSamples);
1862 jgs 119 if (isTagged()) {
1863     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1864     tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1865     }
1866     }
1867    
1868     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1869     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1870     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1871    
1872     //
1873     // Flatten Shape to an array of integers suitable for writing to file
1874     int flatShape[4] = {0,0,0,0};
1875     cout << "\tshape: < ";
1876     for (int dim=0; dim<dataPointRank; dim++) {
1877     flatShape[dim] = dataPointShape[dim];
1878     cout << dataPointShape[dim] << " ";
1879     }
1880     cout << ">" << endl;
1881    
1882     //
1883 jgs 123 // Open archive file
1884 jgs 119 ofstream archiveFile;
1885     archiveFile.open(fileName.data(), ios::out);
1886    
1887     if (!archiveFile.good()) {
1888     throw DataException("archiveData Error: problem opening archive file");
1889     }
1890    
1891 jgs 123 //
1892     // Write common data items to archive file
1893 jgs 119 archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1894     archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1895     archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1896     archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1897     archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1898     archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1899     archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1900     for (int dim = 0; dim < 4; dim++) {
1901     archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1902     }
1903     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1904     archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1905     }
1906     if (isTagged()) {
1907     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1908     archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1909     }
1910     }
1911    
1912     if (!archiveFile.good()) {
1913     throw DataException("archiveData Error: problem writing to archive file");
1914     }
1915    
1916     //
1917 jgs 123 // Archive underlying data values for each Data type
1918     int noValues;
1919 jgs 119 switch (dataType) {
1920     case 0:
1921     // DataEmpty
1922 jgs 123 noValues = 0;
1923     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1924     cout << "\tnoValues: " << noValues << endl;
1925 jgs 119 break;
1926     case 1:
1927     // DataConstant
1928 jgs 123 noValues = m_data->getLength();
1929     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1930     cout << "\tnoValues: " << noValues << endl;
1931     if (m_data->archiveData(archiveFile,noValues)) {
1932     throw DataException("archiveData Error: problem writing data to archive file");
1933     }
1934 jgs 119 break;
1935     case 2:
1936     // DataTagged
1937 jgs 123 noValues = m_data->getLength();
1938     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1939     cout << "\tnoValues: " << noValues << endl;
1940     if (m_data->archiveData(archiveFile,noValues)) {
1941     throw DataException("archiveData Error: problem writing data to archive file");
1942     }
1943 jgs 119 break;
1944     case 3:
1945     // DataExpanded
1946 jgs 123 noValues = m_data->getLength();
1947     archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1948     cout << "\tnoValues: " << noValues << endl;
1949     if (m_data->archiveData(archiveFile,noValues)) {
1950     throw DataException("archiveData Error: problem writing data to archive file");
1951     }
1952 jgs 119 break;
1953     }
1954    
1955 jgs 123 if (!archiveFile.good()) {
1956     throw DataException("archiveData Error: problem writing data to archive file");
1957     }
1958    
1959     //
1960     // Close archive file
1961     archiveFile.close();
1962    
1963     if (!archiveFile.good()) {
1964     throw DataException("archiveData Error: problem closing archive file");
1965     }
1966    
1967 jgs 119 }
1968    
1969     void
1970     Data::extractData(const std::string fileName,
1971     const FunctionSpace& fspace)
1972     {
1973     //
1974     // Can only extract Data to an object which is initially DataEmpty
1975     if (!isEmpty()) {
1976     throw DataException("extractData Error: can only extract to DataEmpty object");
1977     }
1978    
1979     cout << "Extracting Data object from: " << fileName << endl;
1980    
1981     int dataType;
1982     int noSamples;
1983     int noDPPSample;
1984     int functionSpaceType;
1985     int dataPointRank;
1986     int dataPointSize;
1987     int dataLength;
1988     DataArrayView::ShapeType dataPointShape;
1989     int flatShape[4];
1990    
1991     //
1992 jgs 123 // Open the archive file
1993 jgs 119 ifstream archiveFile;
1994     archiveFile.open(fileName.data(), ios::in);
1995    
1996     if (!archiveFile.good()) {
1997     throw DataException("extractData Error: problem opening archive file");
1998     }
1999    
2000 jgs 123 //
2001     // Read common data items from archive file
2002 jgs 119 archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
2003     archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2004     archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2005     archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2006     archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2007     archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2008     archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
2009     for (int dim = 0; dim < 4; dim++) {
2010     archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2011     if (flatShape[dim]>0) {
2012     dataPointShape.push_back(flatShape[dim]);
2013     }
2014     }
2015 woo409 757 vector<int> referenceNumbers(noSamples);
2016 jgs 119 for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2017     archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2018     }
2019 woo409 757 vector<int> tagNumbers(noSamples);
2020 jgs 119 if (dataType==2) {
2021     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2022     archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2023     }
2024     }
2025    
2026     if (!archiveFile.good()) {
2027     throw DataException("extractData Error: problem reading from archive file");
2028     }
2029    
2030 jgs 123 //
2031     // Verify the values just read from the archive file
2032 jgs 119 switch (dataType) {
2033     case 0:
2034     cout << "\tdataType: DataEmpty" << endl;
2035     break;
2036     case 1:
2037     cout << "\tdataType: DataConstant" << endl;
2038     break;
2039     case 2:
2040     cout << "\tdataType: DataTagged" << endl;
2041     break;
2042     case 3:
2043     cout << "\tdataType: DataExpanded" << endl;
2044     break;
2045     default:
2046     throw DataException("extractData Error: undefined dataType read from archive file");
2047     break;
2048     }
2049    
2050     cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2051     cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2052     cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2053     cout << "\tshape: < ";
2054     for (int dim = 0; dim < dataPointRank; dim++) {
2055     cout << dataPointShape[dim] << " ";
2056     }
2057     cout << ">" << endl;
2058    
2059     //
2060     // Verify that supplied FunctionSpace object is compatible with this Data object.
2061     if ( (fspace.getTypeCode()!=functionSpaceType) ||
2062     (fspace.getNumSamples()!=noSamples) ||
2063     (fspace.getNumDPPSample()!=noDPPSample)
2064     ) {
2065     throw DataException("extractData Error: incompatible FunctionSpace");
2066     }
2067     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2068 gross 964 if (referenceNumbers[sampleNo] != fspace.getReferenceIDOfSample(sampleNo)) {
2069 jgs 119 throw DataException("extractData Error: incompatible FunctionSpace");
2070     }
2071     }
2072     if (dataType==2) {
2073     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2074     if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2075     throw DataException("extractData Error: incompatible FunctionSpace");
2076     }
2077     }
2078     }
2079    
2080     //
2081     // Construct a DataVector to hold underlying data values
2082     DataVector dataVec(dataLength);
2083    
2084     //
2085     // Load this DataVector with the appropriate values
2086 jgs 123 int noValues;
2087     archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2088     cout << "\tnoValues: " << noValues << endl;
2089 jgs 119 switch (dataType) {
2090     case 0:
2091     // DataEmpty
2092 jgs 123 if (noValues != 0) {
2093     throw DataException("extractData Error: problem reading data from archive file");
2094     }
2095 jgs 119 break;
2096     case 1:
2097     // DataConstant
2098 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2099     throw DataException("extractData Error: problem reading data from archive file");
2100     }
2101 jgs 119 break;
2102     case 2:
2103     // DataTagged
2104 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2105     throw DataException("extractData Error: problem reading data from archive file");
2106     }
2107 jgs 119 break;
2108     case 3:
2109     // DataExpanded
2110 jgs 123 if (dataVec.extractData(archiveFile,noValues)) {
2111     throw DataException("extractData Error: problem reading data from archive file");
2112     }
2113 jgs 119 break;
2114     }
2115    
2116 jgs 123 if (!archiveFile.good()) {
2117     throw DataException("extractData Error: problem reading from archive file");
2118     }
2119    
2120 jgs 119 //
2121 jgs 123 // Close archive file
2122     archiveFile.close();
2123    
2124     if (!archiveFile.good()) {
2125     throw DataException("extractData Error: problem closing archive file");
2126     }
2127    
2128     //
2129 jgs 119 // Construct an appropriate Data object
2130     DataAbstract* tempData;
2131     switch (dataType) {
2132     case 0:
2133     // DataEmpty
2134     tempData=new DataEmpty();
2135     break;
2136     case 1:
2137     // DataConstant
2138     tempData=new DataConstant(fspace,dataPointShape,dataVec);
2139     break;
2140     case 2:
2141     // DataTagged
2142     tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2143     break;
2144     case 3:
2145     // DataExpanded
2146     tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2147     break;
2148     }
2149     shared_ptr<DataAbstract> temp_data(tempData);
2150     m_data=temp_data;
2151     }
2152    
2153 jgs 94 ostream& escript::operator<<(ostream& o, const Data& data)
2154     {
2155     o << data.toString();
2156     return o;
2157     }
2158 bcumming 782
2159 ksteube 813 Data
2160     escript::C_GeneralTensorProduct(Data& arg_0,
2161     Data& arg_1,
2162     int axis_offset,
2163     int transpose)
2164     {
2165 gross 826 // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2166 ksteube 813 // SM is the product of the last axis_offset entries in arg_0.getShape().
2167    
2168     // Interpolate if necessary and find an appropriate function space
2169 gross 826 Data arg_0_Z, arg_1_Z;
2170 ksteube 813 if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2171     if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2172 gross 826 arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2173     arg_1_Z = Data(arg_1);
2174 ksteube 813 }
2175     else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2176 gross 826 arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2177     arg_0_Z =Data(arg_0);
2178 ksteube 813 }
2179     else {
2180     throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2181     }
2182 gross 826 } else {
2183     arg_0_Z = Data(arg_0);
2184     arg_1_Z = Data(arg_1);
2185 ksteube 813 }
2186     // Get rank and shape of inputs
2187 gross 826 int rank0 = arg_0_Z.getDataPointRank();
2188     int rank1 = arg_1_Z.getDataPointRank();
2189     DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2190     DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2191 ksteube 813
2192     // Prepare for the loops of the product and verify compatibility of shapes
2193     int start0=0, start1=0;
2194     if (transpose == 0) {}
2195     else if (transpose == 1) { start0 = axis_offset; }
2196     else if (transpose == 2) { start1 = rank1-axis_offset; }
2197     else { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2198    
2199     // Adjust the shapes for transpose
2200     DataArrayView::ShapeType tmpShape0;
2201     DataArrayView::ShapeType tmpShape1;
2202     for (int i=0; i<rank0; i++) { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2203     for (int i=0; i<rank1; i++) { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2204    
2205     #if 0
2206     // For debugging: show shape after transpose
2207     char tmp[100];
2208     std::string shapeStr;
2209     shapeStr = "(";
2210     for (int i=0; i<rank0; i++) { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2211     shapeStr += ")";
2212     cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2213     shapeStr = "(";
2214     for (int i=0; i<rank1; i++) { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2215     shapeStr += ")";
2216     cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2217     #endif
2218    
2219     // Prepare for the loops of the product
2220     int SL=1, SM=1, SR=1;
2221     for (int i=0; i<rank0-axis_offset; i++) {
2222     SL *= tmpShape0[i];
2223     }
2224     for (int i=rank0-axis_offset; i<rank0; i++) {
2225     if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2226     throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2227     }
2228     SM *= tmpShape0[i];
2229     }
2230     for (int i=axis_offset; i<rank1; i++) {
2231     SR *= tmpShape1[i];
2232     }
2233    
2234     // Define the shape of the output
2235     DataArrayView::ShapeType shape2;
2236 gross 826 for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2237     for (int i=axis_offset; i<rank1; i++) { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2238 ksteube 813
2239     // Declare output Data object
2240 gross 826 Data res;
2241 ksteube 813
2242 gross 826 if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2243     res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataConstant output
2244     double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2245     double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2246     double *ptr_2 = &((res.getPointDataView().getData())[0]);
2247 ksteube 813 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2248     }
2249 gross 826 else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2250 ksteube 813
2251     // Prepare the DataConstant input
2252 gross 826 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2253 ksteube 813 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2254    
2255     // Borrow DataTagged input from Data object
2256 gross 826 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2257 ksteube 813 if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2258    
2259     // Prepare a DataTagged output 2
2260 gross 826 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataTagged output
2261     res.tag();
2262     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2263 ksteube 813 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2264    
2265     // Prepare offset into DataConstant
2266     int offset_0 = tmp_0->getPointOffset(0,0);
2267 gross 826 double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2268 ksteube 813 // Get the views
2269     DataArrayView view_1 = tmp_1->getDefaultValue();
2270     DataArrayView view_2 = tmp_2->getDefaultValue();
2271     // Get the pointers to the actual data
2272     double *ptr_1 = &((view_1.getData())[0]);
2273     double *ptr_2 = &((view_2.getData())[0]);
2274     // Compute an MVP for the default
2275     matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2276     // Compute an MVP for each tag
2277     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2278     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2279     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2280     tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2281     DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2282     DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2283     double *ptr_1 = &view_1.getData(0);
2284     double *ptr_2 = &view_2.getData(0);
2285     matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2286     }
2287    
2288     }
2289 gross 826 else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2290 ksteube 813
2291 gross 826 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2292     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2293     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2294     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2295 ksteube 813 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2296     if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2297     if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2298     int sampleNo_1,dataPointNo_1;
2299 gross 826 int numSamples_1 = arg_1_Z.getNumSamples();
2300     int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2301 ksteube 813 int offset_0 = tmp_0->getPointOffset(0,0);
2302     #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2303     for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2304     for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2305     int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2306     int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2307 gross 826 double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2308     double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2309     double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2310 ksteube 813 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2311     }
2312     }
2313    
2314     }
2315 gross 826 else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2316 ksteube 813
2317     // Borrow DataTagged input from Data object
2318 gross 826 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2319 ksteube 813 if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2320    
2321     // Prepare the DataConstant input
2322 gross 826 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2323 ksteube 813 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2324    
2325     // Prepare a DataTagged output 2
2326 gross 826 res = Data(0.0, shape2, arg_0_Z.getFunctionSpace()); // DataTagged output
2327     res.tag();
2328     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2329 ksteube 813 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2330    
2331     // Prepare offset into DataConstant
2332     int offset_1 = tmp_1->getPointOffset(0,0);
2333 gross 826 double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2334 ksteube 813 // Get the views
2335     DataArrayView view_0 = tmp_0->getDefaultValue();
2336     DataArrayView view_2 = tmp_2->getDefaultValue();
2337     // Get the pointers to the actual data
2338     double *ptr_0 = &((view_0.getData())[0]);
2339     double *ptr_2 = &((view_2.getData())[0]);
2340     // Compute an MVP for the default
2341     matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2342     // Compute an MVP for each tag
2343     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2344     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2345     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2346     tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2347     DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2348     DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2349     double *ptr_0 = &view_0.getData(0);
2350     double *ptr_2 = &view_2.getData(0);
2351     matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2352     }
2353    
2354     }
2355 gross 826 else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2356 ksteube 813
2357     // Borrow DataTagged input from Data object
2358 gross 826 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2359 ksteube 813 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2360    
2361     // Borrow DataTagged input from Data object
2362 gross 826 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2363 ksteube 813 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2364    
2365     // Prepare a DataTagged output 2
2366 gross 826 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2367     res.tag(); // DataTagged output
2368     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2369 ksteube 813 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2370    
2371     // Get the views
2372     DataArrayView view_0 = tmp_0->getDefaultValue();
2373     DataArrayView view_1 = tmp_1->getDefaultValue();
2374     DataArrayView view_2 = tmp_2->getDefaultValue();
2375     // Get the pointers to the actual data
2376     double *ptr_0 = &((view_0.getData())[0]);
2377     double *ptr_1 = &((view_1.getData())[0]);
2378     double *ptr_2 = &((view_2.getData())[0]);
2379     // Compute an MVP for the default
2380     matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2381     // Merge the tags
2382     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2383     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2384     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2385     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2386     tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2387     }
2388     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2389     tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2390     }
2391     // Compute an MVP for each tag
2392     const