/[escript]/trunk-mpi-branch/escript/src/Data.cpp
ViewVC logotype

Annotation of /trunk-mpi-branch/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1308 - (hide annotations)
Wed Sep 19 06:39:45 2007 UTC (12 years, 1 month ago) by matt
File size: 75618 byte(s)
Removed redundant code.

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