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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1092 - (show annotations)
Fri Apr 13 03:39:49 2007 UTC (12 years, 2 months ago) by gross
File size: 74839 byte(s)
the useless profilinf for data.cpp removed (doDebug=yes should work again)
and a small bug in the gmsh reader fixed.

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26