/[escript]/branches/arrayview_from_1695_trunk/escript/src/Data.cpp
ViewVC logotype

Annotation of /branches/arrayview_from_1695_trunk/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1712 - (hide annotations)
Wed Aug 20 05:04:28 2008 UTC (12 years, 2 months ago) by jfenwick
File size: 77255 byte(s)
Branch commit.

Finished first pass of Data.h - There is still a constructor which takes 
a DataArrayView as a parameter. Apart from that, there are no direct 
references to DataArrayView.

DataTagged has a new constructor for copying just the tags from an 
existing object.
DataTypes:: now has a scalarShape constant (to avoid creating one 
everytime you create a scalar).



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