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

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

Parent Directory Parent Directory | Revision Log Revision Log


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