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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26