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

Diff of /trunk/escript/src/Data.h

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/esys2/escript/src/Data/Data.h revision 97 by jgs, Tue Dec 14 05:39:33 2004 UTC trunk/escript/src/Data.h revision 2459 by jfenwick, Thu Jun 4 06:17:54 2009 UTC
# Line 1  Line 1 
 // $Id$  
 /*=============================================================================  
1    
2   ******************************************************************************  /*******************************************************
3   *                                                                            *  *
4   *       COPYRIGHT ACcESS 2004 -  All Rights Reserved                         *  * Copyright (c) 2003-2008 by University of Queensland
5   *                                                                            *  * Earth Systems Science Computational Center (ESSCC)
6   * This software is the property of ACcESS.  No part of this code             *  * http://www.uq.edu.au/esscc
7   * may be copied in any form or by any means without the expressed written    *  *
8   * consent of ACcESS.  Copying, use or modification of this software          *  * Primary Business: Queensland, Australia
9   * by any unauthorised person is illegal unless that                          *  * Licensed under the Open Software License version 3.0
10   * person has a software license agreement with ACcESS.                       *  * http://www.opensource.org/licenses/osl-3.0.php
11   *                                                                            *  *
12   ******************************************************************************  *******************************************************/
13    
14  ******************************************************************************/  
15    /** \file Data.h */
16    
17  #ifndef DATA_H  #ifndef DATA_H
18  #define DATA_H  #define DATA_H
19    #include "system_dep.h"
20    
21    #include "DataTypes.h"
22    #include "DataAbstract.h"
23    #include "DataAlgorithm.h"
24    #include "FunctionSpace.h"
25    #include "BinaryOp.h"
26    #include "UnaryOp.h"
27    #include "DataException.h"
28    
 #include "escript/Data/DataAbstract.h"  
 #include "escript/Data/DataTagged.h"  
 #include "escript/Data/FunctionSpace.h"  
 #include "escript/Data/BinaryOp.h"  
 #include "escript/Data/UnaryOp.h"  
29    
30  extern "C" {  extern "C" {
31  #include "escript/Data/DataC.h"  #include "DataC.h"
32    //#include <omp.h>
33  }  }
34    
35  #include <iostream>  #include "esysmpi.h"
36  #include <string>  #include <string>
 #include <memory>  
37  #include <algorithm>  #include <algorithm>
38    #include <sstream>
39    
40  #include <boost/shared_ptr.hpp>  #include <boost/shared_ptr.hpp>
41  #include <boost/python/object.hpp>  #include <boost/python/object.hpp>
 #include <boost/python/list.hpp>  
42  #include <boost/python/tuple.hpp>  #include <boost/python/tuple.hpp>
 #include <boost/python/numeric.hpp>  
43    
44  /**  #include "BufferGroup.h"
    \brief  
    Data is essentially a factory class which creates the appropriate Data  
    object for the given construction arguments. It retains control over  
    the object created for the lifetime of the object.  
    The type of Data object referred to may change during the lifetime of  
    the Data object.  
   
    Description:  
    Data is essentially a factory class which creates the appropriate Data  
    object for the given construction arguments. It retains control over  
    the object created for the lifetime of the object.  
    The type of Data object referred to may change during the lifetime of  
    the Data object.  
 */  
45    
46  namespace escript {  namespace escript {
47    
48  //  //
49  // Forward declaration for various implimentations of Data.  // Forward declaration for various implementations of Data.
 class DataEmpty;  
50  class DataConstant;  class DataConstant;
51  class DataTagged;  class DataTagged;
52  class DataExpanded;  class DataExpanded;
53    class DataLazy;
54    
55    /**
56       \brief
57       Data represents a collection of datapoints.
58    
59       Description:
60       Internally, the datapoints are actually stored by a DataAbstract object.
61       The specific instance of DataAbstract used may vary over the lifetime
62       of the Data object.
63       Some methods on this class return references (eg getShape()).
64       These references should not be used after an operation which changes the underlying DataAbstract object.
65       Doing so will lead to invalid memory access.
66       This should not affect any methods exposed via boost::python.
67    */
68  class Data {  class Data {
69    
70    public:    public:
71    
72      // These typedefs allow function names to be cast to pointers
73      // to functions of the appropriate type when calling unaryOp etc.
74    typedef double (*UnaryDFunPtr)(double);    typedef double (*UnaryDFunPtr)(double);
75    typedef double (*BinaryDFunPtr)(double,double);    typedef double (*BinaryDFunPtr)(double,double);
76    
77    
78    /**    /**
79       Constructors.       Constructors.
80    */    */
# Line 80  class Data { Line 84  class Data {
84       Default constructor.       Default constructor.
85       Creates a DataEmpty object.       Creates a DataEmpty object.
86    */    */
87      ESCRIPT_DLL_API
88    Data();    Data();
89    
90    /**    /**
# Line 87  class Data { Line 92  class Data {
92       Copy constructor.       Copy constructor.
93       WARNING: Only performs a shallow copy.       WARNING: Only performs a shallow copy.
94    */    */
95      ESCRIPT_DLL_API
96    Data(const Data& inData);    Data(const Data& inData);
97    
98    /**    /**
# Line 95  class Data { Line 101  class Data {
101       function space of inData the inData are tried to be interpolated to what,       function space of inData the inData are tried to be interpolated to what,
102       otherwise a shallow copy of inData is returned.       otherwise a shallow copy of inData is returned.
103    */    */
104      ESCRIPT_DLL_API
105    Data(const Data& inData,    Data(const Data& inData,
106         const FunctionSpace& what);         const FunctionSpace& what);
107    
108    /**    /**
109       \brief      \brief Copy Data from an existing vector
110       Constructor which copies data from a DataArrayView.    */
111    
112       \param value - Input - Data value for a single point.    ESCRIPT_DLL_API
113       \param what - Input - A description of what this data represents.    Data(const DataTypes::ValueType& value,
114       \param expanded - Input - Flag, if true fill the entire container with           const DataTypes::ShapeType& shape,
115                         the value. Otherwise a more efficient storage                   const FunctionSpace& what=FunctionSpace(),
116                         mechanism will be used.                   bool expanded=false);
   */  
   Data(const DataArrayView& value,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
117    
118    /**    /**
119       \brief       \brief
120       Constructor which creates a Data from a DataArrayView shape.       Constructor which creates a Data with points having the specified shape.
121    
122       \param value - Input - Single value applied to all Data.       \param value - Input - Single value applied to all Data.
123       \param dataPointShape - Input - The shape of each data point.       \param dataPointShape - Input - The shape of each data point.
# Line 123  class Data { Line 126  class Data {
126                         the given value. Otherwise a more efficient storage                         the given value. Otherwise a more efficient storage
127                         mechanism will be used.                         mechanism will be used.
128    */    */
129      ESCRIPT_DLL_API
130    Data(double value,    Data(double value,
131         const DataArrayView::ShapeType& dataPointShape=DataArrayView::ShapeType(),         const DataTypes::ShapeType& dataPointShape=DataTypes::ShapeType(),
132         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
133         bool expanded=false);         bool expanded=false);
134    
# Line 135  class Data { Line 139  class Data {
139       \param inData - Input - Input Data object.       \param inData - Input - Input Data object.
140       \param region - Input - Region to copy.       \param region - Input - Region to copy.
141    */    */
142      ESCRIPT_DLL_API
143    Data(const Data& inData,    Data(const Data& inData,
144         const DataArrayView::RegionType& region);         const DataTypes::RegionType& region);
   
   /**  
      \brief  
      Constructor which will create Tagged data if expanded is false.  
      No attempt is made to ensure the tag keys match the tag keys  
      within the function space.  
   
      \param tagKeys - Input - List of tag values.  
      \param values - Input - List of values, one for each tag.  
      \param defaultValue - Input - A default value, used if tag doesn't exist.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the appropriate values.  
   */  
   Data(const DataTagged::TagListType& tagKeys,  
        const DataTagged::ValueListType& values,  
        const DataArrayView& defaultValue,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
   
   /**  
      \brief  
      Constructor which copies data from a python numarray.  
   
      \param value - Input - Data value for a single point.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the value. Otherwise a more efficient storage  
                        mechanism will be used.  
   */  
   Data(const boost::python::numeric::array& value,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
145    
146    /**    /**
147       \brief       \brief
148       Constructor which copies data from any object that can be converted into       Constructor which copies data from any object that can be treated like a python array/sequence.
      a python numarray.  
149    
150       \param value - Input - Input data.       \param value - Input - Input data.
151       \param what - Input - A description of what this data represents.       \param what - Input - A description of what this data represents.
# Line 182  class Data { Line 153  class Data {
153                         the value. Otherwise a more efficient storage                         the value. Otherwise a more efficient storage
154                         mechanism will be used.                         mechanism will be used.
155    */    */
156      ESCRIPT_DLL_API
157    Data(const boost::python::object& value,    Data(const boost::python::object& value,
158         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
159         bool expanded=false);         bool expanded=false);
# Line 189  class Data { Line 161  class Data {
161    /**    /**
162       \brief       \brief
163       Constructor which creates a DataConstant.       Constructor which creates a DataConstant.
164       Copies data from any object that can be converted       Copies data from any object that can be treated like a python array/sequence.
165       into a numarray. All other parameters are copied from other.       All other parameters are copied from other.
166    
167       \param value - Input - Input data.       \param value - Input - Input data.
168       \param other - Input - contains all other parameters.       \param other - Input - contains all other parameters.
169    */    */
170      ESCRIPT_DLL_API
171    Data(const boost::python::object& value,    Data(const boost::python::object& value,
172         const Data& other);         const Data& other);
173    
# Line 202  class Data { Line 175  class Data {
175       \brief       \brief
176       Constructor which creates a DataConstant of "shape" with constant value.       Constructor which creates a DataConstant of "shape" with constant value.
177    */    */
178    Data(double value,    ESCRIPT_DLL_API
179         const boost::python::tuple& shape=boost::python::make_tuple(),    Data(double value,
180           const boost::python::tuple& shape=boost::python::make_tuple(),
181         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
182         bool expanded=false);         bool expanded=false);
183    
184    
185    
186      /**
187        \brief Create a Data using an existing DataAbstract. Warning: The new object assumes ownership of the pointer!
188        Once you have passed the pointer, do not delete it.
189      */
190      ESCRIPT_DLL_API
191      explicit Data(DataAbstract* underlyingdata);
192    
193      /**
194        \brief Create a Data based on the supplied DataAbstract
195      */
196      ESCRIPT_DLL_API
197      explicit Data(DataAbstract_ptr underlyingdata);
198    
199    /**    /**
200       \brief       \brief
201       Perform a deep copy.       Destructor
202      */
203      ESCRIPT_DLL_API
204      ~Data();
205    
206      /**
207         \brief Make this object a deep copy of "other".
208    */    */
209      ESCRIPT_DLL_API
210    void    void
211    copy(const Data& other);    copy(const Data& other);
212    
213    /**    /**
214         \brief Return a pointer to a deep copy of this object.
215      */
216      ESCRIPT_DLL_API
217      Data
218      copySelf();
219    
220    
221      /**
222         \brief produce a delayed evaluation version of this Data.
223      */
224      ESCRIPT_DLL_API
225      Data
226      delay();
227    
228      /**
229         \brief convert the current data into lazy data.
230      */
231      ESCRIPT_DLL_API
232      void
233      delaySelf();
234    
235    
236      /**
237       Member access methods.       Member access methods.
238    */    */
239    
240    /**    /**
241       \brief       \brief
242         switches on update protection
243    
244      */
245      ESCRIPT_DLL_API
246      void
247      setProtection();
248    
249      /**
250         \brief
251         Returns trueif the data object is protected against update
252    
253      */
254      ESCRIPT_DLL_API
255      bool
256      isProtected() const;
257    
258    
259    /**
260       \brief
261       Return teh value of a data point as a python tuple.
262    */
263      ESCRIPT_DLL_API
264      const boost::python::object
265      getValueOfDataPointAsTuple(int dataPointNo);
266    
267      /**
268         \brief
269         sets the values of a data-point from a python object on this process
270      */
271      ESCRIPT_DLL_API
272      void
273      setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object);
274    
275      /**
276         \brief
277         sets the values of a data-point from a array-like object on this process
278      */
279      ESCRIPT_DLL_API
280      void
281      setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);
282    
283      /**
284         \brief
285         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 Return a data point across all processors as a python tuple.
293      */
294      ESCRIPT_DLL_API
295      const boost::python::object
296      getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo);
297    
298      /**
299         \brief
300         Return the tag number associated with the given data-point.
301    
302      */
303      ESCRIPT_DLL_API
304      int
305      getTagNumber(int dpno);
306    
307      /**
308         \brief
309       Return the C wrapper for the Data object.       Return the C wrapper for the Data object.
310    */    */
311      ESCRIPT_DLL_API
312    escriptDataC    escriptDataC
313    getDataC();    getDataC();
314    
315    
316    
317    /**    /**
318       \brief       \brief
319       Return the C wrapper for the Data object - const version.       Return the C wrapper for the Data object - const version.
320    */    */
321      ESCRIPT_DLL_API
322    escriptDataC    escriptDataC
323    getDataC() const;    getDataC() const;
324    
325    /**    /**
326       \brief      \brief How much space is required to evaulate a sample of the Data.
      Write the data as a string.  
327    */    */
328    inline    ESCRIPT_DLL_API
329    std::string    size_t
330    toString() const    getSampleBufferSize() const;
331    {  
332      return m_data->toString();  
   }  
333    
334    /**    /**
335       \brief       \brief
336       Return the DataArrayView of the point data. This essentially contains       Write the data as a string. For large amounts of data, a summary is printed.
      the shape information for each data point although it also may be used  
      to manipulate the point data.  
337    */    */
338    inline    ESCRIPT_DLL_API
339    const DataArrayView&    std::string
340    getPointDataView() const    toString() const;
   {  
      return m_data->getPointDataView();  
   }  
341    
342    /**    /**
343       \brief       \brief
344       Whatever the current Data type make this into a DataExpanded.       Whatever the current Data type make this into a DataExpanded.
345    */    */
346      ESCRIPT_DLL_API
347    void    void
348    expand();    expand();
349    
# Line 269  class Data { Line 353  class Data {
353       Constant data to be converted to tagged. An attempt to convert       Constant data to be converted to tagged. An attempt to convert
354       Expanded data to tagged will throw an exception.       Expanded data to tagged will throw an exception.
355    */    */
356      ESCRIPT_DLL_API
357    void    void
358    tag();    tag();
359    
360    /**    /**
361        \brief If this data is lazy, then convert it to ready data.
362        What type of ready data depends on the expression. For example, Constant+Tagged==Tagged.
363      */
364      ESCRIPT_DLL_API
365      void
366      resolve();
367    
368    
369      /**
370       \brief Ensures data is ready for write access.
371      This means that the data will be resolved if lazy and will be copied if shared with another Data object.
372      \warning This method should only be called in single threaded sections of code. (It modifies m_data).
373      Do not create any Data objects from this one between calling requireWrite and getSampleDataRW.
374      Doing so might introduce additional sharing.
375      */
376      ESCRIPT_DLL_API
377      void
378      requireWrite();
379    
380      /**
381       \brief       \brief
382       Return true if this Data is expanded.       Return true if this Data is expanded.
383         \note To determine if a sample will contain separate values for each datapoint. Use actsExpanded instead.
384    */    */
385      ESCRIPT_DLL_API
386    bool    bool
387    isExpanded() const;    isExpanded() const;
388    
389    /**    /**
390       \brief       \brief
391         Return true if this Data is expanded or resolves to expanded.
392         That is, if it has a separate value for each datapoint in the sample.
393      */
394      ESCRIPT_DLL_API
395      bool
396      actsExpanded() const;
397      
398    
399      /**
400         \brief
401       Return true if this Data is tagged.       Return true if this Data is tagged.
402    */    */
403      ESCRIPT_DLL_API
404    bool    bool
405    isTagged() const;    isTagged() const;
406    
# Line 290  class Data { Line 408  class Data {
408       \brief       \brief
409       Return true if this Data is constant.       Return true if this Data is constant.
410    */    */
411      ESCRIPT_DLL_API
412    bool    bool
413    isConstant() const;    isConstant() const;
414    
415    /**    /**
416         \brief Return true if this Data is lazy.
417      */
418      ESCRIPT_DLL_API
419      bool
420      isLazy() const;
421    
422      /**
423         \brief Return true if this data is ready.
424      */
425      ESCRIPT_DLL_API
426      bool
427      isReady() const;
428    
429      /**
430       \brief       \brief
431       Return true if this Data is empty.       Return true if this Data holds an instance of DataEmpty. This is _not_ the same as asking if the object
432    contains datapoints.
433    */    */
434      ESCRIPT_DLL_API
435    bool    bool
436    isEmpty() const;    isEmpty() const;
437    
# Line 304  class Data { Line 439  class Data {
439       \brief       \brief
440       Return the function space.       Return the function space.
441    */    */
442      ESCRIPT_DLL_API
443    inline    inline
444    const FunctionSpace&    const FunctionSpace&
445    getFunctionSpace() const    getFunctionSpace() const
# Line 315  class Data { Line 451  class Data {
451       \brief       \brief
452       Return a copy of the function space.       Return a copy of the function space.
453    */    */
454      ESCRIPT_DLL_API
455    const FunctionSpace    const FunctionSpace
456    getCopyOfFunctionSpace() const;    getCopyOfFunctionSpace() const;
457    
# Line 322  class Data { Line 459  class Data {
459       \brief       \brief
460       Return the domain.       Return the domain.
461    */    */
462      ESCRIPT_DLL_API
463    inline    inline
464    const AbstractDomain&  //   const AbstractDomain&
465      const_Domain_ptr
466    getDomain() const    getDomain() const
467    {    {
468       return getFunctionSpace().getDomain();       return getFunctionSpace().getDomain();
469    }    }
470    
471    
472      /**
473         \brief
474         Return the domain.
475         TODO: For internal use only.   This should be removed.
476      */
477      ESCRIPT_DLL_API
478      inline
479    //   const AbstractDomain&
480      Domain_ptr
481      getDomainPython() const
482      {
483         return getFunctionSpace().getDomainPython();
484      }
485    
486    /**    /**
487       \brief       \brief
488       Return a copy of the domain.       Return a copy of the domain.
489    */    */
490      ESCRIPT_DLL_API
491    const AbstractDomain    const AbstractDomain
492    getCopyOfDomain() const;    getCopyOfDomain() const;
493    
# Line 340  class Data { Line 495  class Data {
495       \brief       \brief
496       Return the rank of the point data.       Return the rank of the point data.
497    */    */
498      ESCRIPT_DLL_API
499    inline    inline
500    int    unsigned int
501    getDataPointRank() const    getDataPointRank() const
502    {    {
503      return m_data->getPointDataView().getRank();      return m_data->getRank();
504    }    }
505    
506    /**    /**
507       \brief       \brief
508         Return the number of data points
509      */
510      ESCRIPT_DLL_API
511      inline
512      int
513      getNumDataPoints() const
514      {
515        return getNumSamples() * getNumDataPointsPerSample();
516      }
517      /**
518         \brief
519       Return the number of samples.       Return the number of samples.
520    */    */
521      ESCRIPT_DLL_API
522    inline    inline
523    int    int
524    getNumSamples() const    getNumSamples() const
# Line 362  class Data { Line 530  class Data {
530       \brief       \brief
531       Return the number of data points per sample.       Return the number of data points per sample.
532    */    */
533      ESCRIPT_DLL_API
534    inline    inline
535    int    int
536    getNumDataPointsPerSample() const    getNumDataPointsPerSample() const
# Line 369  class Data { Line 538  class Data {
538      return m_data->getNumDPPSample();      return m_data->getNumDPPSample();
539    }    }
540    
541    
542      /**
543        \brief
544        Return the number of values in the shape for this object.
545      */
546      ESCRIPT_DLL_API
547      int
548      getNoValues() const
549      {
550        return m_data->getNoValues();
551      }
552    
553    
554      /**
555         \brief
556         dumps the object into a netCDF file
557      */
558      ESCRIPT_DLL_API
559      void
560      dump(const std::string fileName) const;
561    
562     /**
563      \brief returns the values of the object as a list of tuples (one for each datapoint).
564    
565      \param scalarastuple If true, scalar data will produce single valued tuples [(1,) (2,) ...]
566    If false, the result is a list of scalars [1, 2, ...]
567     */
568      ESCRIPT_DLL_API
569      const boost::python::object
570      toListOfTuples(bool scalarastuple=false);
571    
572    
573     /**
574        \brief
575        Return the sample data for the given sample no. This is not the
576        preferred interface but is provided for use by C code.
577        The bufferg parameter is only required for LazyData.
578        \param sampleNo - Input - the given sample no.
579        \param buffer - A buffer to to compute (and store) sample data in will be selected from this group.
580        \return pointer to the sample data.
581    */
582      ESCRIPT_DLL_API
583      inline
584      const DataAbstract::ValueType::value_type*
585      getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg=0);
586    
587    
588    /**    /**
589       \brief       \brief
590       Return the sample data for the given sample no. This is not the       Return the sample data for the given sample no. This is not the
591       preferred interface but is provided for use by C code.       preferred interface but is provided for use by C code.
592       \param sampleNo - Input - the given sample no.       \param sampleNo - Input - the given sample no.
593         \return pointer to the sample data.
594    */    */
595      ESCRIPT_DLL_API
596    inline    inline
597    DataAbstract::ValueType::value_type*    DataAbstract::ValueType::value_type*
598    getSampleData(DataAbstract::ValueType::size_type sampleNo)    getSampleDataRW(DataAbstract::ValueType::size_type sampleNo);
599    {  
     return m_data->getSampleData(sampleNo);  
   }  
600    
601    /**    /**
602       \brief       \brief
# Line 388  class Data { Line 604  class Data {
604       access data that isn't tagged an exception will be thrown.       access data that isn't tagged an exception will be thrown.
605       \param tag - Input - the tag key.       \param tag - Input - the tag key.
606    */    */
607      ESCRIPT_DLL_API
608    inline    inline
609    DataAbstract::ValueType::value_type*    DataAbstract::ValueType::value_type*
610    getSampleDataByTag(int tag)    getSampleDataByTag(int tag)
# Line 397  class Data { Line 614  class Data {
614    
615    /**    /**
616       \brief       \brief
617       Return a view into the data for the data point specified.       Return a reference into the DataVector which points to the specified data point.
618       NOTE: Construction of the DataArrayView is a relatively expensive       \param sampleNo - Input -
619       operation.       \param dataPointNo - Input -
620      */
621      ESCRIPT_DLL_API
622      DataTypes::ValueType::const_reference
623      getDataPointRO(int sampleNo, int dataPointNo);
624    
625      /**
626         \brief
627         Return a reference into the DataVector which points to the specified data point.
628       \param sampleNo - Input -       \param sampleNo - Input -
629       \param dataPointNo - Input -       \param dataPointNo - Input -
630    */    */
631      ESCRIPT_DLL_API
632      DataTypes::ValueType::reference
633      getDataPointRW(int sampleNo, int dataPointNo);
634    
635    
636    
637      /**
638         \brief
639         Return the offset for the given sample and point within the sample
640      */
641      ESCRIPT_DLL_API
642    inline    inline
643    DataArrayView    DataTypes::ValueType::size_type
644    getDataPoint(int sampleNo,    getDataOffset(int sampleNo,
645                 int dataPointNo)                 int dataPointNo)
646    {    {
647      return m_data->getDataPoint(sampleNo,dataPointNo);        return m_data->getPointOffset(sampleNo,dataPointNo);
648    }    }
649    
650    /**    /**
651       \brief       \brief
652       Return a reference to the data point shape.       Return a reference to the data point shape.
653    */    */
654    const DataArrayView::ShapeType&    ESCRIPT_DLL_API
655    getDataPointShape() const;    inline
656      const DataTypes::ShapeType&
657      getDataPointShape() const
658      {
659        return m_data->getShape();
660      }
661    
662    /**    /**
663       \brief       \brief
664       Return the data point shape as a tuple of integers.       Return the data point shape as a tuple of integers.
665    */    */
666    boost::python::tuple    ESCRIPT_DLL_API
667      const boost::python::tuple
668    getShapeTuple() const;    getShapeTuple() const;
669    
670    /**    /**
# Line 430  class Data { Line 672  class Data {
672       Return the size of the data point. It is the product of the       Return the size of the data point. It is the product of the
673       data point shape dimensions.       data point shape dimensions.
674    */    */
675      ESCRIPT_DLL_API
676    int    int
677    getDataPointSize() const;    getDataPointSize() const;
678    
# Line 437  class Data { Line 680  class Data {
680       \brief       \brief
681       Return the number of doubles stored for this Data.       Return the number of doubles stored for this Data.
682    */    */
683    DataArrayView::ValueType::size_type    ESCRIPT_DLL_API
684      DataTypes::ValueType::size_type
685    getLength() const;    getLength() const;
686    
687    
688    
689    /**    /**
690       \brief       \brief
691       Assign the given value to the tag. Implicitly converts this       Assign the given value to the tag assocciated with name. Implicitly converts this
692       object to type DataTagged. Throws an exception if this object       object to type DataTagged. Throws an exception if this object
693       cannot be converted to a DataTagged object.       cannot be converted to a DataTagged object or name cannot be mapped onto a tag key.
694       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
695       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
696        ==>*
697    */    */
698      ESCRIPT_DLL_API
699    void    void
700    setTaggedValue(int tagKey,    setTaggedValueByName(std::string name,
701                   const boost::python::object& value);                         const boost::python::object& value);
702    
703    /**    /**
704       \brief       \brief
705       Assign the given value to the tag. Implicitly converts this       Assign the given value to the tag. Implicitly converts this
706       object to type DataTagged. Throws an exception if this object       object to type DataTagged if it is constant.
707       cannot be converted to a DataTagged object.  
708       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
709       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
710       Note: removed for now - this version not needed, and breaks escript.cpp      ==>*
711    */    */
712    /*    ESCRIPT_DLL_API
713    void    void
714    setTaggedValue(int tagKey,    setTaggedValue(int tagKey,
715                   const DataArrayView& value);                   const boost::python::object& value);
716    
717      /**
718         \brief
719         Assign the given value to the tag. Implicitly converts this
720         object to type DataTagged if it is constant.
721    
722         \param tagKey - Input - Integer key.
723         \param pointshape - Input - The shape of the value parameter
724         \param value - Input - Value to associate with given key.
725         \param dataOffset - Input - Offset of the begining of the point within the value parameter
726    */    */
727      ESCRIPT_DLL_API
728      void
729      setTaggedValueFromCPP(int tagKey,
730                const DataTypes::ShapeType& pointshape,
731                            const DataTypes::ValueType& value,
732                int dataOffset=0);
733    
734    
735    
736    /**    /**
737      \brief      \brief
738      Copy other Data object into this Data object where mask is positive.      Copy other Data object into this Data object where mask is positive.
739    */    */
740      ESCRIPT_DLL_API
741    void    void
742    copyWithMask(const Data& other,    copyWithMask(const Data& other,
743                 const Data& mask);                 const Data& mask);
# Line 481  class Data { Line 748  class Data {
748    
749    /**    /**
750       \brief       \brief
751         set all values to zero
752         *
753      */
754      ESCRIPT_DLL_API
755      void
756      setToZero();
757    
758      /**
759         \brief
760       Interpolates this onto the given functionspace and returns       Interpolates this onto the given functionspace and returns
761       the result as a Data object.       the result as a Data object.
762         *
763    */    */
764      ESCRIPT_DLL_API
765    Data    Data
766    interpolate(const FunctionSpace& functionspace) const;    interpolate(const FunctionSpace& functionspace) const;
   
767    /**    /**
768       \brief       \brief
769       Calculates the gradient of the data at the data points of functionspace.       Calculates the gradient of the data at the data points of functionspace.
770       If functionspace is not present the function space of Function(getDomain()) is used.       If functionspace is not present the function space of Function(getDomain()) is used.
771         *
772    */    */
773      ESCRIPT_DLL_API
774    Data    Data
775    gradOn(const FunctionSpace& functionspace) const;    gradOn(const FunctionSpace& functionspace) const;
776    
777      ESCRIPT_DLL_API
778    Data    Data
779    grad() const;    grad() const;
780    
781    /**    /**
782       \brief      \brief
783       Calculate the integral over the function space domain.       Calculate the integral over the function space domain as a python tuple.
784      */
785      ESCRIPT_DLL_API
786      boost::python::object
787      integrateToTuple_const() const;
788    
789    
790      /**
791        \brief
792         Calculate the integral over the function space domain as a python tuple.
793    */    */
794    boost::python::numeric::array    ESCRIPT_DLL_API
795    integrate() const;    boost::python::object
796      integrateToTuple();
797    
798    
799    
800    /**    /**
801       \brief       \brief
802         Returns 1./ Data object
803         *
804      */
805      ESCRIPT_DLL_API
806      Data
807      oneOver() const;
808      /**
809         \brief
810       Return a Data with a 1 for +ive values and a 0 for 0 or -ive values.       Return a Data with a 1 for +ive values and a 0 for 0 or -ive values.
811         *
812    */    */
813      ESCRIPT_DLL_API
814    Data    Data
815    wherePositive() const;    wherePositive() const;
816    
817    /**    /**
818       \brief       \brief
819       Return a Data with a 1 for -ive values and a 0 for +ive or 0 values.       Return a Data with a 1 for -ive values and a 0 for +ive or 0 values.
820         *
821    */    */
822      ESCRIPT_DLL_API
823    Data    Data
824    whereNegative() const;    whereNegative() const;
825    
826    /**    /**
827       \brief       \brief
828       Return a Data with a 1 for +ive or 0 values and a 0 for -ive values.       Return a Data with a 1 for +ive or 0 values and a 0 for -ive values.
829         *
830    */    */
831      ESCRIPT_DLL_API
832    Data    Data
833    whereNonNegative() const;    whereNonNegative() const;
834    
835    /**    /**
836       \brief       \brief
837       Return a Data with a 1 for -ive or 0 values and a 0 for +ive values.       Return a Data with a 1 for -ive or 0 values and a 0 for +ive values.
838         *
839    */    */
840      ESCRIPT_DLL_API
841    Data    Data
842    whereNonPositive() const;    whereNonPositive() const;
843    
844    /**    /**
845       \brief       \brief
846       Return a Data with a 1 for 0 values and a 0 for +ive or -ive values.       Return a Data with a 1 for 0 values and a 0 for +ive or -ive values.
847         *
848    */    */
849      ESCRIPT_DLL_API
850    Data    Data
851    whereZero() const;    whereZero(double tol=0.0) const;
852    
853    /**    /**
854       \brief       \brief
855       Return a Data with a 0 for 0 values and a 1 for +ive or -ive values.       Return a Data with a 0 for 0 values and a 1 for +ive or -ive values.
856         *
857      */
858      ESCRIPT_DLL_API
859      Data
860      whereNonZero(double tol=0.0) const;
861    
862      /**
863         \brief
864         Return the maximum absolute value of this Data object.
865    
866         The method is not const because lazy data needs to be expanded before Lsup can be computed.
867         The _const form can be used when the Data object is const, however this will only work for
868         Data which is not Lazy.
869    
870         For Data which contain no samples (or tagged Data for which no tags in use have a value)
871         zero is returned.
872      */
873      ESCRIPT_DLL_API
874      double
875      Lsup();
876    
877      ESCRIPT_DLL_API
878      double
879      Lsup_const() const;
880    
881    
882      /**
883         \brief
884         Return the maximum value of this Data object.
885    
886         The method is not const because lazy data needs to be expanded before sup can be computed.
887         The _const form can be used when the Data object is const, however this will only work for
888         Data which is not Lazy.
889    
890         For Data which contain no samples (or tagged Data for which no tags in use have a value)
891         a large negative value is returned.
892      */
893      ESCRIPT_DLL_API
894      double
895      sup();
896    
897      ESCRIPT_DLL_API
898      double
899      sup_const() const;
900    
901    
902      /**
903         \brief
904         Return the minimum value of this Data object.
905    
906         The method is not const because lazy data needs to be expanded before inf can be computed.
907         The _const form can be used when the Data object is const, however this will only work for
908         Data which is not Lazy.
909    
910         For Data which contain no samples (or tagged Data for which no tags in use have a value)
911         a large positive value is returned.
912      */
913      ESCRIPT_DLL_API
914      double
915      inf();
916    
917      ESCRIPT_DLL_API
918      double
919      inf_const() const;
920    
921    
922    
923      /**
924         \brief
925         Return the absolute value of each data point of this Data object.
926         *
927      */
928      ESCRIPT_DLL_API
929      Data
930      abs() const;
931    
932      /**
933         \brief
934         Return the maximum value of each data point of this Data object.
935         *
936      */
937      ESCRIPT_DLL_API
938      Data
939      maxval() const;
940    
941      /**
942         \brief
943         Return the minimum value of each data point of this Data object.
944         *
945      */
946      ESCRIPT_DLL_API
947      Data
948      minval() const;
949    
950      /**
951         \brief
952         Return the (sample number, data-point number) of the data point with
953         the minimum value in this Data object.
954      */
955      ESCRIPT_DLL_API
956      const boost::python::tuple
957      minGlobalDataPoint() const;
958    
959      ESCRIPT_DLL_API
960      void
961      calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
962      /**
963         \brief
964         Return the sign of each data point of this Data object.
965         -1 for negative values, zero for zero values, 1 for positive values.
966         *
967      */
968      ESCRIPT_DLL_API
969      Data
970      sign() const;
971    
972      /**
973         \brief
974         Return the symmetric part of a matrix which is half the matrix plus its transpose.
975         *
976      */
977      ESCRIPT_DLL_API
978      Data
979      symmetric() const;
980    
981      /**
982         \brief
983         Return the nonsymmetric part of a matrix which is half the matrix minus its transpose.
984         *
985      */
986      ESCRIPT_DLL_API
987      Data
988      nonsymmetric() const;
989    
990      /**
991         \brief
992         Return the trace of a matrix
993         *
994      */
995      ESCRIPT_DLL_API
996      Data
997      trace(int axis_offset) const;
998    
999      /**
1000         \brief
1001         Transpose each data point of this Data object around the given axis.
1002         *
1003    */    */
1004      ESCRIPT_DLL_API
1005    Data    Data
1006    whereNonZero() const;    transpose(int axis_offset) const;
1007    
1008      /**
1009         \brief
1010         Return the eigenvalues of the symmetric part at each data point of this Data object in increasing values.
1011         Currently this function is restricted to rank 2, square shape, and dimension 3.
1012         *
1013      */
1014      ESCRIPT_DLL_API
1015      Data
1016      eigenvalues() const;
1017    
1018      /**
1019         \brief
1020         Return the eigenvalues and corresponding eigenvcetors of the symmetric part at each data point of this Data object.
1021         the eigenvalues are ordered in increasing size where eigenvalues with relative difference less than
1022         tol are treated as equal. The eigenvectors are orthogonal, normalized and the sclaed such that the
1023         first non-zero entry is positive.
1024         Currently this function is restricted to rank 2, square shape, and dimension 3
1025         *
1026      */
1027      ESCRIPT_DLL_API
1028      const boost::python::tuple
1029      eigenvalues_and_eigenvectors(const double tol=1.e-12) const;
1030    
1031      /**
1032         \brief
1033         swaps the components axis0 and axis1
1034         *
1035      */
1036      ESCRIPT_DLL_API
1037      Data
1038      swapaxes(const int axis0, const int axis1) const;
1039    
1040      /**
1041         \brief
1042         Return the error function erf of each data point of this Data object.
1043         *
1044      */
1045      ESCRIPT_DLL_API
1046      Data
1047      erf() const;
1048    
1049    /**    /**
1050       \brief       \brief
1051       Return the sin of each data point of this Data object.       Return the sin of each data point of this Data object.
1052         *
1053    */    */
1054      ESCRIPT_DLL_API
1055    Data    Data
1056    sin() const;    sin() const;
1057    
1058    /**    /**
1059       \brief       \brief
1060       Return the cos of each data point of this Data object.       Return the cos of each data point of this Data object.
1061         *
1062    */    */
1063      ESCRIPT_DLL_API
1064    Data    Data
1065    cos() const;    cos() const;
1066    
1067    /**    /**
1068       \brief       \brief
1069       Return the tan of each data point of this Data object.       Return the tan of each data point of this Data object.
1070         *
1071    */    */
1072      ESCRIPT_DLL_API
1073    Data    Data
1074    tan() const;    tan() const;
1075    
1076    /**    /**
1077       \brief       \brief
1078       Return the log to base 10 of each data point of this Data object.       Return the asin of each data point of this Data object.
1079         *
1080    */    */
1081      ESCRIPT_DLL_API
1082    Data    Data
1083    log() const;    asin() const;
1084    
1085    /**    /**
1086       \brief       \brief
1087       Return the natural log of each data point of this Data object.       Return the acos of each data point of this Data object.
1088         *
1089    */    */
1090      ESCRIPT_DLL_API
1091    Data    Data
1092    ln() const;    acos() const;
1093    
1094    /**    /**
1095       \brief       \brief
1096       Return the maximum absolute value of this Data object.       Return the atan of each data point of this Data object.
1097         *
1098    */    */
1099    double    ESCRIPT_DLL_API
1100    Lsup() const;    Data
1101      atan() const;
1102    
1103    /**    /**
1104       \brief       \brief
1105       Return the maximum value of this Data object.       Return the sinh of each data point of this Data object.
1106         *
1107    */    */
1108    double    ESCRIPT_DLL_API
1109    sup() const;    Data
1110      sinh() const;
1111    
1112    /**    /**
1113       \brief       \brief
1114       Return the minimum value of this Data object.       Return the cosh of each data point of this Data object.
1115         *
1116    */    */
1117    double    ESCRIPT_DLL_API
1118    inf() const;    Data
1119      cosh() const;
1120    
1121    /**    /**
1122       \brief       \brief
1123       Return the absolute value of each data point of this Data object.       Return the tanh of each data point of this Data object.
1124         *
1125    */    */
1126      ESCRIPT_DLL_API
1127    Data    Data
1128    abs() const;    tanh() const;
1129    
1130    /**    /**
1131       \brief       \brief
1132       Return the maximum value of each data point of this Data object.       Return the asinh of each data point of this Data object.
1133         *
1134    */    */
1135      ESCRIPT_DLL_API
1136    Data    Data
1137    maxval() const;    asinh() const;
1138    
1139    /**    /**
1140       \brief       \brief
1141       Return the minimum value of each data point of this Data object.       Return the acosh of each data point of this Data object.
1142         *
1143    */    */
1144      ESCRIPT_DLL_API
1145    Data    Data
1146    minval() const;    acosh() const;
1147    
1148    /**    /**
1149       \brief       \brief
1150       Return the length of each data point of this Data object.       Return the atanh of each data point of this Data object.
1151       sqrt(sum(A[i,j,k,l]^2))       *
1152    */    */
1153      ESCRIPT_DLL_API
1154    Data    Data
1155    length() const;    atanh() const;
1156    
1157    /**    /**
1158       \brief       \brief
1159       Return the sign of each data point of this Data object.       Return the log to base 10 of each data point of this Data object.
1160       -1 for negative values, zero for zero values, 1 for positive values.       *
1161    */    */
1162      ESCRIPT_DLL_API
1163    Data    Data
1164    sign() const;    log10() const;
1165    
1166    /**    /**
1167      \transpose       \brief
1168      Transpose each data point of this Data object around the given axis.       Return the natural log of each data point of this Data object.
1169         *
1170    */    */
1171      ESCRIPT_DLL_API
1172    Data    Data
1173    transpose(int axis) const;    log() const;
1174    
1175    /**    /**
1176      \trace       \brief
1177      Calculate the trace of each data point of this Data object.       Return the exponential function of each data point of this Data object.
1178      sum(A[i,i,i,i])       *
1179    */    */
1180      ESCRIPT_DLL_API
1181    Data    Data
1182    trace() const;    exp() const;
1183    
1184    /**    /**
1185      \exp       \brief
1186      Return the exponential function of each data point of this Data object.       Return the square root of each data point of this Data object.
1187         *
1188    */    */
1189      ESCRIPT_DLL_API
1190    Data    Data
1191    exp() const;    sqrt() const;
1192    
1193    /**    /**
1194      \sqrt       \brief
1195      Return the square root of each data point of this Data object.       Return the negation of each data point of this Data object.
1196         *
1197    */    */
1198      ESCRIPT_DLL_API
1199    Data    Data
1200    sqrt() const;    neg() const;
1201    
1202      /**
1203         \brief
1204         Return the identity of each data point of this Data object.
1205         Simply returns this object unmodified.
1206         *
1207      */
1208      ESCRIPT_DLL_API
1209      Data
1210      pos() const;
1211    
1212    /**    /**
1213       \brief       \brief
1214       Return the given power of each data point of this Data object.       Return the given power of each data point of this Data object.
1215    
1216         \param right Input - the power to raise the object to.
1217         *
1218    */    */
1219      ESCRIPT_DLL_API
1220    Data    Data
1221    powD(const Data& right) const;    powD(const Data& right) const;
1222    
1223      /**
1224         \brief
1225         Return the given power of each data point of this boost python object.
1226    
1227         \param right Input - the power to raise the object to.
1228         *
1229       */
1230      ESCRIPT_DLL_API
1231    Data    Data
1232    powO(const boost::python::object& right) const;    powO(const boost::python::object& right) const;
1233    
1234    /**    /**
1235      \brief       \brief
1236      Return the negation of each data point of this Data object.       Return the given power of each data point of this boost python object.
1237    */  
1238         \param left Input - the bases
1239         *
1240       */
1241    
1242      ESCRIPT_DLL_API
1243    Data    Data
1244    neg() const;    rpowO(const boost::python::object& left) const;
1245    
1246    /**    /**
1247      \brief       \brief
1248      Return the identity of each data point of this Data object.       writes the object to a file in the DX file format
     Simply returns this object unmodified.  
1249    */    */
1250    Data    ESCRIPT_DLL_API
1251    pos() const;    void
1252      saveDX(std::string fileName) const;
1253    
1254      /**
1255         \brief
1256         writes the object to a file in the VTK file format
1257      */
1258      ESCRIPT_DLL_API
1259      void
1260      saveVTK(std::string fileName) const;
1261    
1262    /**    /**
1263       \brief       \brief
1264       Overloaded operator +=       Overloaded operator +=
1265       \param right - Input - The right hand side.       \param right - Input - The right hand side.
1266         *
1267    */    */
1268      ESCRIPT_DLL_API
1269    Data& operator+=(const Data& right);    Data& operator+=(const Data& right);
1270      ESCRIPT_DLL_API
1271    Data& operator+=(const boost::python::object& right);    Data& operator+=(const boost::python::object& right);
1272    
1273      ESCRIPT_DLL_API
1274      Data& operator=(const Data& other);
1275    
1276    /**    /**
1277       \brief       \brief
1278       Overloaded operator -=       Overloaded operator -=
1279       \param right - Input - The right hand side.       \param right - Input - The right hand side.
1280         *
1281    */    */
1282      ESCRIPT_DLL_API
1283    Data& operator-=(const Data& right);    Data& operator-=(const Data& right);
1284      ESCRIPT_DLL_API
1285    Data& operator-=(const boost::python::object& right);    Data& operator-=(const boost::python::object& right);
1286    
1287   /**   /**
1288       \brief       \brief
1289       Overloaded operator *=       Overloaded operator *=
1290       \param right - Input - The right hand side.       \param right - Input - The right hand side.
1291         *
1292    */    */
1293      ESCRIPT_DLL_API
1294    Data& operator*=(const Data& right);    Data& operator*=(const Data& right);
1295      ESCRIPT_DLL_API
1296    Data& operator*=(const boost::python::object& right);    Data& operator*=(const boost::python::object& right);
1297    
1298   /**   /**
1299       \brief       \brief
1300       Overloaded operator /=       Overloaded operator /=
1301       \param right - Input - The right hand side.       \param right - Input - The right hand side.
1302         *
1303    */    */
1304      ESCRIPT_DLL_API
1305    Data& operator/=(const Data& right);    Data& operator/=(const Data& right);
1306      ESCRIPT_DLL_API
1307    Data& operator/=(const boost::python::object& right);    Data& operator/=(const boost::python::object& right);
1308    
1309    /**    /**
1310       \brief       \brief
1311       Returns true if this can be interpolated to functionspace.       Returns true if this can be interpolated to functionspace.
1312    */    */
1313      ESCRIPT_DLL_API
1314    bool    bool
1315    probeInterpolation(const FunctionSpace& functionspace) const;    probeInterpolation(const FunctionSpace& functionspace) const;
1316    
# Line 748  class Data { Line 1329  class Data {
1329       \param key - Input - python slice tuple specifying       \param key - Input - python slice tuple specifying
1330       slice to return.       slice to return.
1331    */    */
1332      ESCRIPT_DLL_API
1333    Data    Data
1334    getItem(const boost::python::object& key) const;    getItem(const boost::python::object& key) const;
1335    
# Line 755  class Data { Line 1337  class Data {
1337       \brief       \brief
1338       Copies slice from value into this Data object.       Copies slice from value into this Data object.
1339    
      \description  
1340       Implements the [] set operator in python.       Implements the [] set operator in python.
1341       Calls setSlice.       Calls setSlice.
1342    
# Line 763  class Data { Line 1344  class Data {
1344       slice to copy from value.       slice to copy from value.
1345       \param value - Input - Data object to copy from.       \param value - Input - Data object to copy from.
1346    */    */
1347      ESCRIPT_DLL_API
1348    void    void
1349    setItemD(const boost::python::object& key,    setItemD(const boost::python::object& key,
1350             const Data& value);             const Data& value);
1351    
1352      ESCRIPT_DLL_API
1353    void    void
1354    setItemO(const boost::python::object& key,    setItemO(const boost::python::object& key,
1355             const boost::python::object& value);             const boost::python::object& value);
# Line 779  class Data { Line 1362  class Data {
1362       this Data object.       this Data object.
1363    */    */
1364    template <class UnaryFunction>    template <class UnaryFunction>
1365      ESCRIPT_DLL_API
1366    inline    inline
1367    void    void
1368    unaryOp(UnaryFunction operation);    unaryOp2(UnaryFunction operation);
1369    
1370    /**    /**
1371       \brief       \brief
1372       Return a Data object containing the specified slice of       Return a Data object containing the specified slice of
1373       this Data object.       this Data object.
1374       \param region - Input - Region to copy.       \param region - Input - Region to copy.
1375         *
1376    */    */
1377      ESCRIPT_DLL_API
1378    Data    Data
1379    getSlice(const DataArrayView::RegionType& region) const;    getSlice(const DataTypes::RegionType& region) const;
1380    
1381    /**    /**
1382       \brief       \brief
# Line 798  class Data { Line 1384  class Data {
1384       Data object.       Data object.
1385       \param value - Input - Data to copy from.       \param value - Input - Data to copy from.
1386       \param region - Input - Region to copy.       \param region - Input - Region to copy.
1387         *
1388    */    */
1389      ESCRIPT_DLL_API
1390    void    void
1391    setSlice(const Data& value,    setSlice(const Data& value,
1392             const DataArrayView::RegionType& region);             const DataTypes::RegionType& region);
1393    
1394      /**
1395         \brief
1396         print the data values to stdout. Used for debugging
1397      */
1398      ESCRIPT_DLL_API
1399      void
1400            print(void);
1401    
1402      /**
1403         \brief
1404         return the MPI rank number of the local data
1405                     MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()
1406                     is returned
1407      */
1408      ESCRIPT_DLL_API
1409            int
1410            get_MPIRank(void) const;
1411    
1412      /**
1413         \brief
1414         return the MPI rank number of the local data
1415                     MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()
1416                     is returned
1417      */
1418      ESCRIPT_DLL_API
1419            int
1420            get_MPISize(void) const;
1421    
1422      /**
1423         \brief
1424         return the MPI rank number of the local data
1425                     MPI_COMM_WORLD is assumed and returned.
1426      */
1427      ESCRIPT_DLL_API
1428            MPI_Comm
1429            get_MPIComm(void) const;
1430    
1431      /**
1432         \brief
1433         return the object produced by the factory, which is a DataConstant or DataExpanded
1434        TODO Ownership of this object should be explained in doco.
1435      */
1436      ESCRIPT_DLL_API
1437            DataAbstract*
1438            borrowData(void) const;
1439    
1440      ESCRIPT_DLL_API
1441            DataAbstract_ptr
1442            borrowDataPtr(void) const;
1443    
1444      ESCRIPT_DLL_API
1445            DataReady_ptr
1446            borrowReadyPtr(void) const;
1447    
1448    
1449    
1450      /**
1451         \brief
1452         Return a pointer to the beginning of the datapoint at the specified offset.
1453         TODO Eventually these should be inlined.
1454         \param i - position(offset) in the underlying datastructure
1455      */
1456    
1457      ESCRIPT_DLL_API
1458            DataTypes::ValueType::const_reference
1459            getDataAtOffsetRO(DataTypes::ValueType::size_type i);
1460    
1461    
1462      ESCRIPT_DLL_API
1463            DataTypes::ValueType::reference
1464            getDataAtOffsetRW(DataTypes::ValueType::size_type i);
1465    
1466    
1467    
1468    /**
1469       \brief Create a buffer for use by getSample
1470       Allocates a DataVector large enough for DataLazy::resolveSample to operate on for the current Data.
1471       Do not use this buffer for other Data instances (unless you are sure they will be the same size).
1472      
1473       In multi-threaded sections, this needs to be called on each thread.
1474    
1475       \return A BufferGroup* if Data is lazy, NULL otherwise.
1476       \warning This pointer must be deallocated using freeSampleBuffer to avoid cross library memory issues.
1477    */
1478      ESCRIPT_DLL_API
1479      BufferGroup*
1480      allocSampleBuffer() const;
1481    
1482    /**
1483       \brief Free a buffer allocated with allocSampleBuffer.
1484       \param buffer Input - pointer to the buffer to deallocate.
1485    */
1486    ESCRIPT_DLL_API void freeSampleBuffer(BufferGroup* buffer);
1487    
1488   protected:   protected:
1489    
1490   private:   private:
1491    
1492      double
1493      LsupWorker() const;
1494    
1495      double
1496      supWorker() const;
1497    
1498      double
1499      infWorker() const;
1500    
1501      boost::python::object
1502      integrateWorker() const;
1503    
1504    /**    /**
1505       \brief       \brief
1506       Check *this and the right operand are compatible. Throws       Check *this and the right operand are compatible. Throws
# Line 823  class Data { Line 1517  class Data {
1517    /**    /**
1518       \brief       \brief
1519       Perform the specified reduction algorithm on every element of every data point in       Perform the specified reduction algorithm on every element of every data point in
1520       this Data object and return the single double value result.       this Data object according to the given function and return the single value result.
1521    */    */
1522    template <class UnaryFunction>    template <class BinaryFunction>
1523    inline    inline
1524    double    double
1525    algorithm(UnaryFunction operation) const;    algorithm(BinaryFunction operation,
1526                double initial_value) const;
1527    
1528    /**    /**
1529       \brief       \brief
1530       Perform the given binary operation on all of the data's elements.       Reduce each data-point in this Data object using the given operation. Return a Data
1531       The underlying type of the right hand side (right) determines the final       object with the same number of data-points, but with each data-point containing only
1532       type of *this after the operation. For example if the right hand side       one value - the result of the reduction operation on the corresponding data-point in
1533       is expanded *this will be expanded if necessary.       this Data object
      RHS is a Data object.  
1534    */    */
1535    template <class BinaryFunction>    template <class BinaryFunction>
1536    inline    inline
1537    void    Data
1538    binaryOp(const Data& right,    dp_algorithm(BinaryFunction operation,
1539             BinaryFunction operation);                 double initial_value) const;
1540    
1541    /**    /**
1542       \brief       \brief
1543       Perform the given binary operation on all of the data's elements.       Perform the given binary operation on all of the data's elements.
1544       RHS is a boost::python object.       The underlying type of the right hand side (right) determines the final
1545         type of *this after the operation. For example if the right hand side
1546         is expanded *this will be expanded if necessary.
1547         RHS is a Data object.
1548    */    */
1549    template <class BinaryFunction>    template <class BinaryFunction>
1550    inline    inline
1551    void    void
1552    binaryOp(const boost::python::object& right,    binaryOp(const Data& right,
1553             BinaryFunction operation);             BinaryFunction operation);
1554    
1555    /**    /**
# Line 875  class Data { Line 1572  class Data {
1572       \brief       \brief
1573       Construct a Data object of the appropriate type.       Construct a Data object of the appropriate type.
1574    */    */
1575    template <class IValueType>  
1576    void    void
1577    initialise(const IValueType& value,    initialise(const DataTypes::ValueType& value,
1578             const DataTypes::ShapeType& shape,
1579               const FunctionSpace& what,               const FunctionSpace& what,
1580               bool expanded);               bool expanded);
1581    
   /**  
      \brief  
      Reshape the data point if the data point is currently rank 0.  
      Will throw an exception if the data points are not rank 0.  
      The original data point value is used for all values of the new  
      data point.  
   */  
1582    void    void
1583    reshapeDataPoint(const DataArrayView::ShapeType& shape);    initialise(const WrappedArray& value,
1584                     const FunctionSpace& what,
1585                     bool expanded);
1586    
1587      //
1588      // flag to protect the data object against any update
1589      bool m_protected;
1590      mutable bool m_shared;
1591      bool m_lazy;
1592    
1593    //    //
1594    // pointer to the actual data object    // pointer to the actual data object
1595    boost::shared_ptr<DataAbstract> m_data;  //   boost::shared_ptr<DataAbstract> m_data;
1596      DataAbstract_ptr m_data;
1597    
1598  };  // If possible please use getReadyPtr instead.
1599    // But see warning below.
1600      const DataReady*
1601      getReady() const;
1602    
1603  template <class IValueType>    DataReady*
1604  void    getReady();
1605  Data::initialise(const IValueType& value,  
1606                   const FunctionSpace& what,  
1607                   bool expanded)  // Be wary of using this for local operations since it (temporarily) increases reference count.
1608  {  // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1609    // getReady() instead
1610      DataReady_ptr
1611      getReadyPtr();
1612    
1613      const_DataReady_ptr
1614      getReadyPtr() const;
1615    
1616    
1617      /**
1618       \brief Update the Data's shared flag
1619       This indicates that the DataAbstract used by this object is now shared (or no longer shared).
1620       For internal use only.
1621      */
1622      void updateShareStatus(bool nowshared) const
1623      {
1624        m_shared=nowshared;     // m_shared is mutable
1625      }
1626    
1627      // In the isShared() method below:
1628      // A problem would occur if m_data (the address pointed to) were being modified
1629      // while the call m_data->is_shared is being executed.
1630      //
1631      // Q: So why do I think this code can be thread safe/correct?
1632      // A: We need to make some assumptions.
1633      // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1634      // 2. We assume that no constructions or assignments which will share previously unshared
1635      //    will occur while this call is executing. This is consistent with the way Data:: and C are written.
1636    //    //
1637    // Construct a Data object of the appropriate type.    // This means that the only transition we need to consider, is when a previously shared object is
1638    // Construct the object first as there seems to be a bug which causes    // not shared anymore. ie. the other objects have been destroyed or a deep copy has been made.
1639    // undefined behaviour if an exception is thrown during construction    // In those cases the m_shared flag changes to false after m_data has completed changing.
1640    // within the shared_ptr constructor.    // For any threads executing before the flag switches they will assume the object is still shared.
1641    if (expanded) {    bool isShared() const
1642      DataAbstract* temp=new DataExpanded(value,what);    {
1643      boost::shared_ptr<DataAbstract> temp_data(temp);      return m_shared;
1644      m_data=temp_data;  /*  if (m_shared) return true;
1645    } else {      if (m_data->isShared())        
1646      DataAbstract* temp=new DataConstant(value,what);      {                  
1647      boost::shared_ptr<DataAbstract> temp_data(temp);          updateShareStatus(true);
1648      m_data=temp_data;          return true;
1649        }
1650        return false;*/
1651    }    }
1652    
1653      void forceResolve()
1654      {
1655        if (isLazy())
1656        {
1657            #ifdef _OPENMP
1658            if (omp_in_parallel())
1659            {   // Yes this is throwing an exception out of an omp thread which is forbidden.
1660            throw DataException("Please do not call forceResolve() in a parallel region.");
1661            }
1662            #endif
1663            resolve();
1664        }
1665      }
1666    
1667      /**
1668      \brief if another object is sharing out member data make a copy to work with instead.
1669      This code should only be called from single threaded sections of code.
1670      */
1671      void exclusiveWrite()
1672      {
1673    #ifdef _OPENMP
1674      if (omp_in_parallel())
1675      {
1676    // *((int*)0)=17;
1677        throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1678      }
1679    #endif
1680        forceResolve();
1681        if (isShared())
1682        {
1683            DataAbstract* t=m_data->deepCopy();
1684            set_m_data(DataAbstract_ptr(t));
1685        }
1686      }
1687    
1688      /**
1689      \brief checks if caller can have exclusive write to the object
1690      */
1691      void checkExclusiveWrite()
1692      {
1693        if  (isLazy() || isShared())
1694        {
1695            throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
1696        }
1697      }
1698    
1699      /**
1700      \brief Modify the data abstract hosted by this Data object
1701      For internal use only.
1702      Passing a pointer to null is permitted (do this in the destructor)
1703      \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1704      */
1705      void set_m_data(DataAbstract_ptr p);
1706    
1707      friend class DataAbstract;        // To allow calls to updateShareStatus
1708    
1709    };
1710    
1711    }   // end namespace escript
1712    
1713    
1714    // No, this is not supposed to be at the top of the file
1715    // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1716    // so that I can dynamic cast between them below.
1717    #include "DataReady.h"
1718    #include "DataLazy.h"
1719    
1720    namespace escript
1721    {
1722    
1723    inline
1724    const DataReady*
1725    Data::getReady() const
1726    {
1727       const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1728       EsysAssert((dr!=0), "Error - casting to DataReady.");
1729       return dr;
1730    }
1731    
1732    inline
1733    DataReady*
1734    Data::getReady()
1735    {
1736       DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1737       EsysAssert((dr!=0), "Error - casting to DataReady.");
1738       return dr;
1739    }
1740    
1741    // Be wary of using this for local operations since it (temporarily) increases reference count.
1742    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1743    // getReady() instead
1744    inline
1745    DataReady_ptr
1746    Data::getReadyPtr()
1747    {
1748       DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1749       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1750       return dr;
1751  }  }
1752    
1753    
1754    inline
1755    const_DataReady_ptr
1756    Data::getReadyPtr() const
1757    {
1758       const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1759       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1760       return dr;
1761    }
1762    
1763    inline
1764    DataAbstract::ValueType::value_type*
1765    Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1766    {
1767    //    if (isLazy())
1768    //    {
1769    //  resolve();
1770    //    }
1771    //    exclusiveWrite();
1772       if (isLazy())
1773       {
1774        throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1775       }
1776       return getReady()->getSampleDataRW(sampleNo);
1777    }
1778    
1779    // inline
1780    // const DataAbstract::ValueType::value_type*
1781    // Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, DataTypes::ValueType* buffer)
1782    // {
1783    //    DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1784    //    if (l!=0)
1785    //    {
1786    //  size_t offset=0;
1787    //  if (buffer==NULL)
1788    //  {
1789    //      throw DataException("Error, attempt to getSampleDataRO for lazy Data with buffer==NULL");
1790    //  }
1791    //  const DataTypes::ValueType* res=l->resolveSample(*buffer,0,sampleNo,offset);
1792    //  return &((*res)[offset]);
1793    //    }
1794    //    return getReady()->getSampleDataRO(sampleNo);
1795    // }
1796    
1797    inline
1798    const DataAbstract::ValueType::value_type*
1799    Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg)
1800    {
1801       DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1802       if (l!=0)
1803       {
1804        size_t offset=0;
1805        if (bufferg==NULL)
1806        {
1807            throw DataException("Error, attempt to getSampleDataRO for lazy Data with buffer==NULL");
1808        }
1809        const DataTypes::ValueType* res=l->resolveSample(*bufferg,sampleNo,offset);
1810        return &((*res)[offset]);
1811       }
1812       return getReady()->getSampleDataRO(sampleNo);
1813    }
1814    
1815    
1816    
1817    /**
1818       Modify a filename for MPI parallel output to multiple files
1819    */
1820    char *Escript_MPI_appendRankToFileName(const char *, int, int);
1821    
1822  /**  /**
1823     Binary Data object operators.     Binary Data object operators.
1824  */  */
1825    inline double rpow(double x,double y)
1826    {
1827        return pow(y,x);
1828    }
1829    
1830  /**  /**
1831    \brief    \brief
1832    Operator+    Operator+
1833    Takes two Data objects.    Takes two Data objects.
1834  */  */
1835  Data operator+(const Data& left, const Data& right);  ESCRIPT_DLL_API Data operator+(const Data& left, const Data& right);
1836    
1837  /**  /**
1838    \brief    \brief
1839    Operator-    Operator-
1840    Takes two Data objects.    Takes two Data objects.
1841  */  */
1842  Data operator-(const Data& left, const Data& right);  ESCRIPT_DLL_API Data operator-(const Data& left, const Data& right);
1843    
1844  /**  /**
1845    \brief    \brief
1846    Operator*    Operator*
1847    Takes two Data objects.    Takes two Data objects.
1848  */  */
1849  Data operator*(const Data& left, const Data& right);  ESCRIPT_DLL_API Data operator*(const Data& left, const Data& right);
1850    
1851  /**  /**
1852    \brief    \brief
1853    Operator/    Operator/
1854    Takes two Data objects.    Takes two Data objects.
1855  */  */
1856  Data operator/(const Data& left, const Data& right);  ESCRIPT_DLL_API Data operator/(const Data& left, const Data& right);
1857    
1858  /**  /**
1859    \brief    \brief
# Line 957  Data operator/(const Data& left, const D Line 1861  Data operator/(const Data& left, const D
1861    Takes LHS Data object and RHS python::object.    Takes LHS Data object and RHS python::object.
1862    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1863  */  */
1864  Data operator+(const Data& left, const boost::python::object& right);  ESCRIPT_DLL_API Data operator+(const Data& left, const boost::python::object& right);
1865    
1866  /**  /**
1867    \brief    \brief
# Line 965  Data operator+(const Data& left, const b Line 1869  Data operator+(const Data& left, const b
1869    Takes LHS Data object and RHS python::object.    Takes LHS Data object and RHS python::object.
1870    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1871  */  */
1872  Data operator-(const Data& left, const boost::python::object& right);  ESCRIPT_DLL_API Data operator-(const Data& left, const boost::python::object& right);
1873    
1874  /**  /**
1875    \brief    \brief
# Line 973  Data operator-(const Data& left, const b Line 1877  Data operator-(const Data& left, const b
1877    Takes LHS Data object and RHS python::object.    Takes LHS Data object and RHS python::object.
1878    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1879  */  */
1880  Data operator*(const Data& left, const boost::python::object& right);  ESCRIPT_DLL_API Data operator*(const Data& left, const boost::python::object& right);
1881    
1882  /**  /**
1883    \brief    \brief
# Line 981  Data operator*(const Data& left, const b Line 1885  Data operator*(const Data& left, const b
1885    Takes LHS Data object and RHS python::object.    Takes LHS Data object and RHS python::object.
1886    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1887  */  */
1888  Data operator/(const Data& left, const boost::python::object& right);  ESCRIPT_DLL_API Data operator/(const Data& left, const boost::python::object& right);
1889    
1890  /**  /**
1891    \brief    \brief
# Line 989  Data operator/(const Data& left, const b Line 1893  Data operator/(const Data& left, const b
1893    Takes LHS python::object and RHS Data object.    Takes LHS python::object and RHS Data object.
1894    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1895  */  */
1896  Data operator+(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator+(const boost::python::object& left, const Data& right);
1897    
1898  /**  /**
1899    \brief    \brief
# Line 997  Data operator+(const boost::python::obje Line 1901  Data operator+(const boost::python::obje
1901    Takes LHS python::object and RHS Data object.    Takes LHS python::object and RHS Data object.
1902    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1903  */  */
1904  Data operator-(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator-(const boost::python::object& left, const Data& right);
1905    
1906  /**  /**
1907    \brief    \brief
# Line 1005  Data operator-(const boost::python::obje Line 1909  Data operator-(const boost::python::obje
1909    Takes LHS python::object and RHS Data object.    Takes LHS python::object and RHS Data object.
1910    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1911  */  */
1912  Data operator*(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator*(const boost::python::object& left, const Data& right);
1913    
1914  /**  /**
1915    \brief    \brief
# Line 1013  Data operator*(const boost::python::obje Line 1917  Data operator*(const boost::python::obje
1917    Takes LHS python::object and RHS Data object.    Takes LHS python::object and RHS Data object.
1918    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1919  */  */
1920  Data operator/(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator/(const boost::python::object& left, const Data& right);
1921    
1922    
1923    
1924  /**  /**
1925    \brief    \brief
1926    Output operator    Output operator
1927  */  */
1928  std::ostream& operator<<(std::ostream& o, const Data& data);  ESCRIPT_DLL_API std::ostream& operator<<(std::ostream& o, const Data& data);
1929    
1930  /**  /**
1931    \brief    \brief
1932    Return true if operands are equivalent, else return false.    Compute a tensor product of two Data objects
1933    NB: this operator does very little at this point, and isn't to    \param arg0 - Input - Data object
1934    be relied on. Requires further implementation.    \param arg1 - Input - Data object
1935      \param axis_offset - Input - axis offset
1936      \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1937  */  */
1938  //bool operator==(const Data& left, const Data& right);  ESCRIPT_DLL_API
1939    Data
1940    C_GeneralTensorProduct(Data& arg0,
1941                         Data& arg1,
1942                         int axis_offset=0,
1943                         int transpose=0);
1944    
1945  /**  /**
1946    \brief    \brief
# Line 1042  Data::binaryOp(const Data& right, Line 1955  Data::binaryOp(const Data& right,
1955  {  {
1956     //     //
1957     // if this has a rank of zero promote it to the rank of the RHS     // if this has a rank of zero promote it to the rank of the RHS
1958     if (getPointDataView().getRank()==0 && right.getPointDataView().getRank()!=0) {     if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1959       reshapeDataPoint(right.getPointDataView().getShape());       throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
1960       }
1961    
1962       if (isLazy() || right.isLazy())
1963       {
1964         throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
1965     }     }
1966     //     //
1967     // initially make the temporary a shallow copy     // initially make the temporary a shallow copy
1968     Data tempRight(right);     Data tempRight(right);
1969    
1970     if (getFunctionSpace()!=right.getFunctionSpace()) {     if (getFunctionSpace()!=right.getFunctionSpace()) {
1971       if (right.probeInterpolation(getFunctionSpace())) {       if (right.probeInterpolation(getFunctionSpace())) {
1972         //         //
1973         // an interpolation is required so create a new Data         // an interpolation is required so create a new Data
1974         tempRight=Data(right,this->getFunctionSpace());         tempRight=Data(right,this->getFunctionSpace());
1975       } else if (probeInterpolation(right.getFunctionSpace())) {       } else if (probeInterpolation(right.getFunctionSpace())) {
1976         //         //
1977         // interpolate onto the RHS function space         // interpolate onto the RHS function space
1978         Data tempLeft(*this,right.getFunctionSpace());         Data tempLeft(*this,right.getFunctionSpace());
1979         m_data=tempLeft.m_data;  //        m_data=tempLeft.m_data;
1980           set_m_data(tempLeft.m_data);
1981       }       }
1982     }     }
1983     operandCheck(tempRight);     operandCheck(tempRight);
# Line 1073  Data::binaryOp(const Data& right, Line 1993  Data::binaryOp(const Data& right,
1993       // of any data type       // of any data type
1994       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
1995       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
1996       escript::binaryOp(*leftC,*(tempRight.m_data.get()),operation);       escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
1997     } else if (isTagged()) {     } else if (isTagged()) {
1998       //       //
1999       // Tagged data is operated on serially, the right hand side can be       // Tagged data is operated on serially, the right hand side can be
# Line 1099  Data::binaryOp(const Data& right, Line 2019  Data::binaryOp(const Data& right,
2019    
2020  /**  /**
2021    \brief    \brief
   Perform the given binary operation with this and right as operands.  
   Right is a boost::python object.  
 */  
 template <class BinaryFunction>  
 inline  
 void  
 Data::binaryOp(const boost::python::object& right,  
                BinaryFunction operation)  
 {  
    DataArray temp(right);  
    //  
    // if this has a rank of zero promote it to the rank of the RHS.  
    if (getPointDataView().getRank()==0 && temp.getView().getRank()!=0) {  
       reshapeDataPoint(temp.getView().getShape());  
    }  
    //  
    // Always allow scalar values for the RHS but check other shapes  
    if (temp.getView().getRank()!=0) {  
      if (!getPointDataView().checkShape(temp.getView().getShape())) {  
        throw DataException(getPointDataView().createShapeErrorMessage(  
                   "Error - RHS shape doesn't match LHS shape.",temp.getView().getShape()));  
      }  
    }  
    if (isExpanded()) {  
      DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());  
      EsysAssert((leftC!=0),"Programming error - casting to DataExpanded.");  
      escript::binaryOp(*leftC,temp.getView(),operation);  
    } else if (isTagged()) {  
      DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());  
      EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");  
      escript::binaryOp(*leftC,temp.getView(),operation);  
    } else if (isConstant()) {  
      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());  
      EsysAssert((leftC!=0),"Programming error - casting to DataConstant.");  
      escript::binaryOp(*leftC,temp.getView(),operation);  
    }  
 }  
   
 /**  
   \brief  
   Perform the given unary operation on other and return the result.  
   Given operation is performed on each element of each data point, thus  
   argument object is a rank n Data object, and returned object is a rank n  
   Data object.  
   Calls Data::unaryOp.  
 */  
 template <class UnaryFunction>  
 inline  
 Data  
 unaryOp(const Data& other,  
         UnaryFunction operation)  
 {  
   Data result;  
   result.copy(other);  
   result.unaryOp(operation);  
   return result;  
 }  
   
 /**  
   \brief  
   Perform the given unary operation on this.  
   Given operation is performed on each element of each data point.  
   Calls escript::unaryOp.  
 */  
 template <class UnaryFunction>  
 inline  
 void  
 Data::unaryOp(UnaryFunction operation)  
 {  
   if (isExpanded()) {  
     DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());  
     EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");  
     escript::unaryOp(*leftC,operation);  
   } else if (isTagged()) {  
     DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());  
     EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");  
     escript::unaryOp(*leftC,operation);  
   } else if (isConstant()) {  
     DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());  
     EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");  
     escript::unaryOp(*leftC,operation);  
   }  
 }  
   
 /**  
   \brief  
2022    Perform the given Data object reduction algorithm on this and return the result.    Perform the given Data object reduction algorithm on this and return the result.
2023    Given operation combines each element of each data point, thus argument    Given operation combines each element of each data point, thus argument
2024    object (*this) is a rank n Data object, and returned object is a scalar.    object (*this) is a rank n Data object, and returned object is a scalar.
2025    Calls escript::algorithm.    Calls escript::algorithm.
2026  */  */
2027  template <class UnaryFunction>  template <class BinaryFunction>
2028  inline  inline
2029  double  double
2030  Data::algorithm(UnaryFunction operation) const  Data::algorithm(BinaryFunction operation, double initial_value) const
2031  {  {
2032    if (isExpanded()) {    if (isExpanded()) {
2033      DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2034      EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");      EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2035      return escript::algorithm(*leftC,operation);      return escript::algorithm(*leftC,operation,initial_value);
2036    } else if (isTagged()) {    } else if (isTagged()) {
2037      DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
2038      EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");      EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");
2039      return escript::algorithm(*leftC,operation);      return escript::algorithm(*leftC,operation,initial_value);
2040    } else if (isConstant()) {    } else if (isConstant()) {
2041      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2042      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2043      return escript::algorithm(*leftC,operation);      return escript::algorithm(*leftC,operation,initial_value);
2044      } else if (isEmpty()) {
2045        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2046      } else if (isLazy()) {
2047        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2048      } else {
2049        throw DataException("Error - Data encapsulates an unknown type.");
2050    }    }
   return 0;  
2051  }  }
2052    
2053  /**  /**
2054    \brief    \brief
2055    Perform the given data point reduction algorithm on data and return the result.    Perform the given data point reduction algorithm on data and return the result.
2056    Given operation combines each element within each data point into a scalar,    Given operation combines each element within each data point into a scalar,
2057    thus argument object is a rank n Data object, and returned object is a    thus argument object is a rank n Data object, and returned object is a
2058    rank 0 Data object.    rank 0 Data object.
2059    Calls escript::dp_algorithm.    Calls escript::dp_algorithm.
2060  */  */
2061  template <class UnaryFunction>  template <class BinaryFunction>
2062  inline  inline
2063  Data  Data
2064  dp_algorithm(const Data& data,  Data::dp_algorithm(BinaryFunction operation, double initial_value) const
              UnaryFunction operation)  
2065  {  {
2066    Data result(0,DataArrayView::ShapeType(),data.getFunctionSpace(),data.isExpanded());    if (isEmpty()) {
2067    if (data.isExpanded()) {      throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2068      DataExpanded* dataE=dynamic_cast<DataExpanded*>(data.m_data.get());    }
2069      else if (isExpanded()) {
2070        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2071        DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2072      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2073      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2074      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2075      escript::dp_algorithm(*dataE,*resultE,operation);      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2076    } else if (data.isTagged()) {      return result;
2077      DataTagged* dataT=dynamic_cast<DataTagged*>(data.m_data.get());    }
2078      DataTagged* resultT=dynamic_cast<DataTagged*>(result.m_data.get());    else if (isTagged()) {
2079        DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());
2080      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2081      EsysAssert((resultT!=0), "Programming error - casting result to DataTagged.");      DataTypes::ValueType defval(1);
2082      escript::dp_algorithm(*dataT,*resultT,operation);      defval[0]=0;
2083    } else if (data.isConstant()) {      DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2084      DataConstant* dataC=dynamic_cast<DataConstant*>(data.m_data.get());      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2085        return Data(resultT);   // note: the Data object now owns the resultT pointer
2086      }
2087      else if (isConstant()) {
2088        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2089        DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2090      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2091      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2092      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2093      escript::dp_algorithm(*dataC,*resultC,operation);      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2094        return result;
2095      } else if (isLazy()) {
2096        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2097      } else {
2098        throw DataException("Error - Data encapsulates an unknown type.");
2099      }
2100    }
2101    
2102    /**
2103      \brief
2104      Compute a tensor operation with two Data objects
2105      \param arg0 - Input - Data object
2106      \param arg1 - Input - Data object
2107      \param operation - Input - Binary op functor
2108    */
2109    template <typename BinaryFunction>
2110    inline
2111    Data
2112    C_TensorBinaryOperation(Data const &arg_0,
2113                            Data const &arg_1,
2114                            BinaryFunction operation)
2115    {
2116      if (arg_0.isEmpty() || arg_1.isEmpty())
2117      {
2118         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2119      }
2120      if (arg_0.isLazy() || arg_1.isLazy())
2121      {
2122         throw DataException("Error - Operations not permitted on lazy data.");
2123      }
2124      // Interpolate if necessary and find an appropriate function space
2125      Data arg_0_Z, arg_1_Z;
2126      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2127        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2128          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2129          arg_1_Z = Data(arg_1);
2130        }
2131        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2132          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2133          arg_0_Z =Data(arg_0);
2134        }
2135        else {
2136          throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible function spaces.");
2137        }
2138      } else {
2139          arg_0_Z = Data(arg_0);
2140          arg_1_Z = Data(arg_1);
2141      }
2142      // Get rank and shape of inputs
2143      int rank0 = arg_0_Z.getDataPointRank();
2144      int rank1 = arg_1_Z.getDataPointRank();
2145      DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2146      DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2147      int size0 = arg_0_Z.getDataPointSize();
2148      int size1 = arg_1_Z.getDataPointSize();
2149      // Declare output Data object
2150      Data res;
2151    
2152      if (shape0 == shape1) {
2153        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2154          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2155          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2156          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2157          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2158    
2159          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2160        }
2161        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2162    
2163          // Prepare the DataConstant input
2164          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2165    
2166          // Borrow DataTagged input from Data object
2167          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2168    
2169          // Prepare a DataTagged output 2
2170          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2171          res.tag();
2172          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2173    
2174          // Prepare offset into DataConstant
2175          int offset_0 = tmp_0->getPointOffset(0,0);
2176          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2177    
2178          // Get the pointers to the actual data
2179          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2180          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2181    
2182          // Compute a result for the default
2183          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2184          // Compute a result for each tag
2185          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2186          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2187          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2188            tmp_2->addTag(i->first);
2189            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2190            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2191    
2192            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2193          }
2194    
2195        }
2196        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2197          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2198          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2199          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2200          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2201    
2202          int sampleNo_1,dataPointNo_1;
2203          int numSamples_1 = arg_1_Z.getNumSamples();
2204          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2205          int offset_0 = tmp_0->getPointOffset(0,0);
2206          res.requireWrite();
2207          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2208          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2209            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2210              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2211              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2212              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2213              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2214              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2215              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2216            }
2217          }
2218    
2219        }
2220        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2221          // Borrow DataTagged input from Data object
2222          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2223    
2224          // Prepare the DataConstant input
2225          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2226    
2227          // Prepare a DataTagged output 2
2228          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2229          res.tag();
2230          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2231    
2232          // Prepare offset into DataConstant
2233          int offset_1 = tmp_1->getPointOffset(0,0);
2234    
2235          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2236          // Get the pointers to the actual data
2237          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2238          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2239          // Compute a result for the default
2240          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2241          // Compute a result for each tag
2242          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2243          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2244          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2245            tmp_2->addTag(i->first);
2246            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2247            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2248            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2249          }
2250    
2251        }
2252        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2253          // Borrow DataTagged input from Data object
2254          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2255    
2256          // Borrow DataTagged input from Data object
2257          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2258    
2259          // Prepare a DataTagged output 2
2260          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2261          res.tag();        // DataTagged output
2262          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2263    
2264          // Get the pointers to the actual data
2265          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2266          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2267          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2268    
2269          // Compute a result for the default
2270          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2271          // Merge the tags
2272          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2273          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2274          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2275          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2276            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2277          }
2278          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2279            tmp_2->addTag(i->first);
2280          }
2281          // Compute a result for each tag
2282          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2283          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2284    
2285            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2286            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2287            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2288    
2289            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2290          }
2291    
2292        }
2293        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2294          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2295          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2296          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(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_0,dataPointNo_0;
2301          int numSamples_0 = arg_0_Z.getNumSamples();
2302          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2303          res.requireWrite();
2304          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2305          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2306            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2307            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2308            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2309              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2310              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2311              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2312              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2313              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2314            }
2315          }
2316    
2317        }
2318        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2319          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2320          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2321          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2322          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2323    
2324          int sampleNo_0,dataPointNo_0;
2325          int numSamples_0 = arg_0_Z.getNumSamples();
2326          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2327          int offset_1 = tmp_1->getPointOffset(0,0);
2328          res.requireWrite();
2329          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2330          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2331            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2332              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2333              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2334    
2335              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2336              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2337              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2338    
2339    
2340              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2341            }
2342          }
2343    
2344        }
2345        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2346          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2347          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2348          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2349          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2350          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2351    
2352          int sampleNo_0,dataPointNo_0;
2353          int numSamples_0 = arg_0_Z.getNumSamples();
2354          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2355          res.requireWrite();
2356          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2357          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2358            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2359            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2360            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2361              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2362              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2363              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2364              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2365              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2366            }
2367          }
2368    
2369        }
2370        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2371          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2372          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2373          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2374          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2375          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2376    
2377          int sampleNo_0,dataPointNo_0;
2378          int numSamples_0 = arg_0_Z.getNumSamples();
2379          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2380          res.requireWrite();
2381          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2382          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2383            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2384              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2385              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2386              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2387              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2388              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2389              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2390              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2391            }
2392          }
2393    
2394        }
2395        else {
2396          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2397        }
2398    
2399      } else if (0 == rank0) {
2400        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2401          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
2402          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2403          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2404          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2405          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2406        }
2407        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2408    
2409          // Prepare the DataConstant input
2410          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2411    
2412          // Borrow DataTagged input from Data object
2413          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2414    
2415          // Prepare a DataTagged output 2
2416          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataTagged output
2417          res.tag();
2418          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2419    
2420          // Prepare offset into DataConstant
2421          int offset_0 = tmp_0->getPointOffset(0,0);
2422          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2423    
2424          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2425          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2426    
2427          // Compute a result for the default
2428          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2429          // Compute a result for each tag
2430          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2431          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2432          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2433            tmp_2->addTag(i->first);
2434            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2435            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2436            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2437          }
2438    
2439        }
2440        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2441    
2442          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2443          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2444          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2445          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2446    
2447          int sampleNo_1,dataPointNo_1;
2448          int numSamples_1 = arg_1_Z.getNumSamples();
2449          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2450          int offset_0 = tmp_0->getPointOffset(0,0);
2451          res.requireWrite();
2452          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2453          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2454            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2455              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2456              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2457              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2458              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2459              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2460              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2461    
2462            }
2463          }
2464    
2465        }
2466        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2467    
2468          // Borrow DataTagged input from Data object
2469          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2470    
2471          // Prepare the DataConstant input
2472          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2473    
2474          // Prepare a DataTagged output 2
2475          res = Data(0.0, shape1, arg_0_Z.getFunctionSpace());      // DataTagged output
2476          res.tag();
2477          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2478    
2479          // Prepare offset into DataConstant
2480          int offset_1 = tmp_1->getPointOffset(0,0);
2481          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2482    
2483          // Get the pointers to the actual data
2484          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2485          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2486    
2487    
2488          // Compute a result for the default
2489          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2490          // Compute a result for each tag
2491          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2492          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2493          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2494            tmp_2->addTag(i->first);
2495            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2496            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2497    
2498            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2499          }
2500    
2501        }
2502        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2503    
2504          // Borrow DataTagged input from Data object
2505          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2506    
2507          // Borrow DataTagged input from Data object
2508          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2509    
2510          // Prepare a DataTagged output 2
2511          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
2512          res.tag();        // DataTagged output
2513          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2514    
2515          // Get the pointers to the actual data
2516          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2517          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2518          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2519    
2520          // Compute a result for the default
2521          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2522          // Merge the tags
2523          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2524          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2525          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2526          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2527            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2528          }
2529          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2530            tmp_2->addTag(i->first);
2531          }
2532          // Compute a result for each tag
2533          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2534          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2535            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2536            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2537            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2538    
2539            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2540          }
2541    
2542        }
2543        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2544    
2545          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2546          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2547          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2548          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2549          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2550    
2551          int sampleNo_0,dataPointNo_0;
2552          int numSamples_0 = arg_0_Z.getNumSamples();
2553          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2554          res.requireWrite();
2555          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2556          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2557            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2558            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2559            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2560              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2561              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2562              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2563              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2564              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2565            }
2566          }
2567    
2568        }
2569        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2570          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2571          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2572          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2573          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2574    
2575          int sampleNo_0,dataPointNo_0;
2576          int numSamples_0 = arg_0_Z.getNumSamples();
2577          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2578          int offset_1 = tmp_1->getPointOffset(0,0);
2579          res.requireWrite();
2580          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2581          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2582            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2583              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2584              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2585              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2586              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2587              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2588              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2589            }
2590          }
2591    
2592    
2593        }
2594        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2595    
2596          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2597          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2598          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2599          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2600          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2601    
2602          int sampleNo_0,dataPointNo_0;
2603          int numSamples_0 = arg_0_Z.getNumSamples();
2604          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2605          res.requireWrite();
2606          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2607          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2608            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2609            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2610            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2611              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2612              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2613              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2614              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2615              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2616            }
2617          }
2618    
2619        }
2620        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2621    
2622          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2623          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2624          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2625          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2626          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2627    
2628          int sampleNo_0,dataPointNo_0;
2629          int numSamples_0 = arg_0_Z.getNumSamples();
2630          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2631          res.requireWrite();
2632          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2633          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2634            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2635              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2636              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2637              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2638              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2639              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2640              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2641              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2642            }
2643          }
2644    
2645        }
2646        else {
2647          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2648        }
2649    
2650      } else if (0 == rank1) {
2651        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2652          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2653          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2654          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2655          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2656          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2657        }
2658        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2659    
2660          // Prepare the DataConstant input
2661          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2662    
2663          // Borrow DataTagged input from Data object
2664          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2665    
2666          // Prepare a DataTagged output 2
2667          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2668          res.tag();
2669          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2670    
2671          // Prepare offset into DataConstant
2672          int offset_0 = tmp_0->getPointOffset(0,0);
2673          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2674    
2675          //Get the pointers to the actual data
2676          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2677          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2678    
2679          // Compute a result for the default
2680          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2681          // Compute a result for each tag
2682          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2683          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2684          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2685            tmp_2->addTag(i->first);
2686            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2687            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2688            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2689          }
2690        }
2691        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2692    
2693          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2694          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2695          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2696          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2697    
2698          int sampleNo_1,dataPointNo_1;
2699          int numSamples_1 = arg_1_Z.getNumSamples();
2700          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2701          int offset_0 = tmp_0->getPointOffset(0,0);
2702          res.requireWrite();
2703          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2704          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2705            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2706              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2707              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2708              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2709              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2710              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2711              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2712            }
2713          }
2714    
2715        }
2716        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2717    
2718          // Borrow DataTagged input from Data object
2719          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2720    
2721          // Prepare the DataConstant input
2722          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2723    
2724          // Prepare a DataTagged output 2
2725          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2726          res.tag();
2727          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2728    
2729          // Prepare offset into DataConstant
2730          int offset_1 = tmp_1->getPointOffset(0,0);
2731          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2732          // Get the pointers to the actual data
2733          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2734          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2735          // Compute a result for the default
2736          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2737          // Compute a result for each tag
2738          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2739          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2740          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2741            tmp_2->addTag(i->first);
2742            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2743            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2744            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2745          }
2746    
2747        }
2748        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2749    
2750          // Borrow DataTagged input from Data object
2751          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2752    
2753          // Borrow DataTagged input from Data object
2754          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2755    
2756          // Prepare a DataTagged output 2
2757          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2758          res.tag();        // DataTagged output
2759          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2760    
2761          // Get the pointers to the actual data
2762          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2763          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2764          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2765    
2766          // Compute a result for the default
2767          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2768          // Merge the tags
2769          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2770          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2771          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2772          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2773            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2774          }
2775          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2776            tmp_2->addTag(i->first);
2777          }
2778          // Compute a result for each tag
2779          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2780          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2781            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2782            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2783            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2784            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2785          }
2786    
2787        }
2788        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2789    
2790          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2791          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2792          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2793          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2794          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2795    
2796          int sampleNo_0,dataPointNo_0;
2797          int numSamples_0 = arg_0_Z.getNumSamples();
2798          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2799          res.requireWrite();
2800          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2801          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2802            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2803            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2804            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2805              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2806              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2807              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2808              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2809              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2810            }
2811          }
2812    
2813        }
2814        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2815          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2816          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2817          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2818          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2819    
2820          int sampleNo_0,dataPointNo_0;
2821          int numSamples_0 = arg_0_Z.getNumSamples();
2822          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2823          int offset_1 = tmp_1->getPointOffset(0,0);
2824          res.requireWrite();
2825          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2826          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2827            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2828              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2829              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2830              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2831              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2832              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2833              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2834            }
2835          }
2836    
2837    
2838        }
2839        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2840    
2841          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2842          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2843          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2844          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2845          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2846    
2847          int sampleNo_0,dataPointNo_0;
2848          int numSamples_0 = arg_0_Z.getNumSamples();
2849          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2850          res.requireWrite();
2851          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2852          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2853            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2854            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2855            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2856              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2857              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2858              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2859              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2860              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2861            }
2862          }
2863    
2864        }
2865        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2866    
2867          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2868          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2869          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2870          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2871          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2872    
2873          int sampleNo_0,dataPointNo_0;
2874          int numSamples_0 = arg_0_Z.getNumSamples();
2875          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2876          res.requireWrite();
2877          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2878          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2879            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2880              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2881              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2882              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2883              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2884              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2885              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2886              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2887            }
2888          }
2889    
2890        }
2891        else {
2892          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2893        }
2894    
2895      } else {
2896        throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible shapes");
2897      }
2898    
2899      return res;
2900    }
2901    
2902    template <typename UnaryFunction>
2903    Data
2904    C_TensorUnaryOperation(Data const &arg_0,
2905                           UnaryFunction operation)
2906    {
2907      if (arg_0.isEmpty())  // do this before we attempt to interpolate
2908      {
2909         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2910      }
2911      if (arg_0.isLazy())
2912      {
2913         throw DataException("Error - Operations not permitted on lazy data.");
2914    }    }
2915    return result;    // Interpolate if necessary and find an appropriate function space
2916      Data arg_0_Z = Data(arg_0);
2917    
2918      // Get rank and shape of inputs
2919      const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2920      int size0 = arg_0_Z.getDataPointSize();
2921    
2922      // Declare output Data object
2923      Data res;
2924    
2925      if (arg_0_Z.isConstant()) {
2926        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
2927        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2928        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2929        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2930      }
2931      else if (arg_0_Z.isTagged()) {
2932    
2933        // Borrow DataTagged input from Data object
2934        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2935    
2936        // Prepare a DataTagged output 2
2937        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());   // DataTagged output
2938        res.tag();
2939        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2940    
2941        // Get the pointers to the actual data
2942        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2943        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2944        // Compute a result for the default
2945        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2946        // Compute a result for each tag
2947        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2948        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2949        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2950          tmp_2->addTag(i->first);
2951          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2952          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2953          tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2954        }
2955    
2956      }
2957      else if (arg_0_Z.isExpanded()) {
2958    
2959        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
2960        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2961        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2962    
2963        int sampleNo_0,dataPointNo_0;
2964        int numSamples_0 = arg_0_Z.getNumSamples();
2965        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2966        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2967        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2968          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2969            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2970            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2971            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2972            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2973            tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2974          }
2975        }
2976      }
2977      else {
2978        throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
2979      }
2980    
2981      return res;
2982  }  }
2983    
2984  }  }

Legend:
Removed from v.97  
changed lines
  Added in v.2459

  ViewVC Help
Powered by ViewVC 1.1.26