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

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

Parent Directory Parent Directory | Revision Log Revision Log


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