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

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

Parent Directory Parent Directory | Revision Log Revision Log


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