/[escript]/branches/schroedinger/escript/src/Data.h
ViewVC logotype

Annotation of /branches/schroedinger/escript/src/Data.h

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26