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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1137 - (show annotations)
Thu May 10 08:11:31 2007 UTC (12 years, 1 month ago) by gross
File size: 75417 byte(s)
This version passes the tests on windows except for 

   * vtk
   * netCDF

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26