/[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 2635 by jfenwick, Thu Aug 27 04:54:41 2009 UTC
# Line 1  Line 1 
 // $Id$  
 /*=============================================================================  
1    
2   ******************************************************************************  /*******************************************************
3   *                                                                            *  *
4   *       COPYRIGHT ACcESS 2004 -  All Rights Reserved                         *  * Copyright (c) 2003-2009 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);
145    
146    /**    /**
147       \brief       \brief
148       Constructor which will create Tagged data if expanded is false.       Constructor which copies data from any object that can be treated like a python array/sequence.
      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);  
   
   /**  
      \brief  
      Constructor which copies data from any object that can be converted into  
      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 true, if 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 the 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       \brief
417       Return true if this Data is empty.       Return true if this Data is constant or resolves to constant.
418         That is, if there is only one value stored, regardless of how many
419         datapoints there are per sample.
420      */
421      ESCRIPT_DLL_API
422      bool
423      actsConstant() const;
424    
425      /**
426         \brief Return true if this Data is lazy.
427    */    */
428      ESCRIPT_DLL_API
429      bool
430      isLazy() const;
431    
432      /**
433         \brief Return true if this data is ready.
434      */
435      ESCRIPT_DLL_API
436      bool
437      isReady() const;
438    
439      /**
440         \brief
441         Return true if this Data holds an instance of DataEmpty. This is _not_ the same as asking if the object
442    contains datapoints.
443      */
444      ESCRIPT_DLL_API
445    bool    bool
446    isEmpty() const;    isEmpty() const;
447    
# Line 304  class Data { Line 449  class Data {
449       \brief       \brief
450       Return the function space.       Return the function space.
451    */    */
452      ESCRIPT_DLL_API
453    inline    inline
454    const FunctionSpace&    const FunctionSpace&
455    getFunctionSpace() const    getFunctionSpace() const
# Line 315  class Data { Line 461  class Data {
461       \brief       \brief
462       Return a copy of the function space.       Return a copy of the function space.
463    */    */
464      ESCRIPT_DLL_API
465    const FunctionSpace    const FunctionSpace
466    getCopyOfFunctionSpace() const;    getCopyOfFunctionSpace() const;
467    
# Line 322  class Data { Line 469  class Data {
469       \brief       \brief
470       Return the domain.       Return the domain.
471    */    */
472      ESCRIPT_DLL_API
473    inline    inline
474    const AbstractDomain&  //   const AbstractDomain&
475      const_Domain_ptr
476    getDomain() const    getDomain() const
477    {    {
478       return getFunctionSpace().getDomain();       return getFunctionSpace().getDomain();
479    }    }
480    
481    
482      /**
483         \brief
484         Return the domain.
485         TODO: For internal use only.   This should be removed.
486      */
487      ESCRIPT_DLL_API
488      inline
489    //   const AbstractDomain&
490      Domain_ptr
491      getDomainPython() const
492      {
493         return getFunctionSpace().getDomainPython();
494      }
495    
496    /**    /**
497       \brief       \brief
498       Return a copy of the domain.       Return a copy of the domain.
499    */    */
500      ESCRIPT_DLL_API
501    const AbstractDomain    const AbstractDomain
502    getCopyOfDomain() const;    getCopyOfDomain() const;
503    
# Line 340  class Data { Line 505  class Data {
505       \brief       \brief
506       Return the rank of the point data.       Return the rank of the point data.
507    */    */
508      ESCRIPT_DLL_API
509    inline    inline
510    int    unsigned int
511    getDataPointRank() const    getDataPointRank() const
512    {    {
513      return m_data->getPointDataView().getRank();      return m_data->getRank();
514    }    }
515    
516    /**    /**
517       \brief       \brief
518         Return the number of data points
519      */
520      ESCRIPT_DLL_API
521      inline
522      int
523      getNumDataPoints() const
524      {
525        return getNumSamples() * getNumDataPointsPerSample();
526      }
527      /**
528         \brief
529       Return the number of samples.       Return the number of samples.
530    */    */
531      ESCRIPT_DLL_API
532    inline    inline
533    int    int
534    getNumSamples() const    getNumSamples() const
# Line 362  class Data { Line 540  class Data {
540       \brief       \brief
541       Return the number of data points per sample.       Return the number of data points per sample.
542    */    */
543      ESCRIPT_DLL_API
544    inline    inline
545    int    int
546    getNumDataPointsPerSample() const    getNumDataPointsPerSample() const
# Line 369  class Data { Line 548  class Data {
548      return m_data->getNumDPPSample();      return m_data->getNumDPPSample();
549    }    }
550    
551    
552      /**
553        \brief
554        Return the number of values in the shape for this object.
555      */
556      ESCRIPT_DLL_API
557      int
558      getNoValues() const
559      {
560        return m_data->getNoValues();
561      }
562    
563    
564      /**
565         \brief
566         dumps the object into a netCDF file
567      */
568      ESCRIPT_DLL_API
569      void
570      dump(const std::string fileName) const;
571    
572     /**
573      \brief returns the values of the object as a list of tuples (one for each datapoint).
574    
575      \param scalarastuple If true, scalar data will produce single valued tuples [(1,) (2,) ...]
576    If false, the result is a list of scalars [1, 2, ...]
577     */
578      ESCRIPT_DLL_API
579      const boost::python::object
580      toListOfTuples(bool scalarastuple=true);
581    
582    
583     /**
584        \brief
585        Return the sample data for the given sample no. This is not the
586        preferred interface but is provided for use by C code.
587        The bufferg parameter is only required for LazyData.
588        \param sampleNo - Input - the given sample no.
589        \param bufferg - A buffer to compute (and store) sample data in will be selected from this group.
590        \return pointer to the sample data.
591    */
592      ESCRIPT_DLL_API
593      inline
594      const DataAbstract::ValueType::value_type*
595      getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg=0);
596    
597    
598    /**    /**
599       \brief       \brief
600       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
601       preferred interface but is provided for use by C code.       preferred interface but is provided for use by C code.
602       \param sampleNo - Input - the given sample no.       \param sampleNo - Input - the given sample no.
603         \return pointer to the sample data.
604    */    */
605      ESCRIPT_DLL_API
606    inline    inline
607    DataAbstract::ValueType::value_type*    DataAbstract::ValueType::value_type*
608    getSampleData(DataAbstract::ValueType::size_type sampleNo)    getSampleDataRW(DataAbstract::ValueType::size_type sampleNo);
609    {  
     return m_data->getSampleData(sampleNo);  
   }  
610    
611    /**    /**
612       \brief       \brief
# Line 388  class Data { Line 614  class Data {
614       access data that isn't tagged an exception will be thrown.       access data that isn't tagged an exception will be thrown.
615       \param tag - Input - the tag key.       \param tag - Input - the tag key.
616    */    */
617      ESCRIPT_DLL_API
618    inline    inline
619    DataAbstract::ValueType::value_type*    DataAbstract::ValueType::value_type*
620    getSampleDataByTag(int tag)    getSampleDataByTag(int tag)
# Line 397  class Data { Line 624  class Data {
624    
625    /**    /**
626       \brief       \brief
627       Return a view into the data for the data point specified.       Return a reference into the DataVector which points to the specified data point.
      NOTE: Construction of the DataArrayView is a relatively expensive  
      operation.  
628       \param sampleNo - Input -       \param sampleNo - Input -
629       \param dataPointNo - Input -       \param dataPointNo - Input -
630    */    */
631      ESCRIPT_DLL_API
632      DataTypes::ValueType::const_reference
633      getDataPointRO(int sampleNo, int dataPointNo);
634    
635      /**
636         \brief
637         Return a reference into the DataVector which points to the specified data point.
638         \param sampleNo - Input -
639         \param dataPointNo - Input -
640      */
641      ESCRIPT_DLL_API
642      DataTypes::ValueType::reference
643      getDataPointRW(int sampleNo, int dataPointNo);
644    
645    
646    
647      /**
648         \brief
649         Return the offset for the given sample and point within the sample
650      */
651      ESCRIPT_DLL_API
652    inline    inline
653    DataArrayView    DataTypes::ValueType::size_type
654    getDataPoint(int sampleNo,    getDataOffset(int sampleNo,
655                 int dataPointNo)                 int dataPointNo)
656    {    {
657      return m_data->getDataPoint(sampleNo,dataPointNo);        return m_data->getPointOffset(sampleNo,dataPointNo);
658    }    }
659    
660    /**    /**
661       \brief       \brief
662       Return a reference to the data point shape.       Return a reference to the data point shape.
663    */    */
664    const DataArrayView::ShapeType&    ESCRIPT_DLL_API
665    getDataPointShape() const;    inline
666      const DataTypes::ShapeType&
667      getDataPointShape() const
668      {
669        return m_data->getShape();
670      }
671    
672    /**    /**
673       \brief       \brief
674       Return the data point shape as a tuple of integers.       Return the data point shape as a tuple of integers.
675    */    */
676    boost::python::tuple    ESCRIPT_DLL_API
677      const boost::python::tuple
678    getShapeTuple() const;    getShapeTuple() const;
679    
680    /**    /**
# Line 430  class Data { Line 682  class Data {
682       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
683       data point shape dimensions.       data point shape dimensions.
684    */    */
685      ESCRIPT_DLL_API
686    int    int
687    getDataPointSize() const;    getDataPointSize() const;
688    
# Line 437  class Data { Line 690  class Data {
690       \brief       \brief
691       Return the number of doubles stored for this Data.       Return the number of doubles stored for this Data.
692    */    */
693    DataArrayView::ValueType::size_type    ESCRIPT_DLL_API
694      DataTypes::ValueType::size_type
695    getLength() const;    getLength() const;
696    
697    
698    
699    /**    /**
700       \brief       \brief
701       Assign the given value to the tag. Implicitly converts this       Assign the given value to the tag assocciated with name. Implicitly converts this
702       object to type DataTagged. Throws an exception if this object       object to type DataTagged. Throws an exception if this object
703       cannot be converted to a DataTagged object.       cannot be converted to a DataTagged object or name cannot be mapped onto a tag key.
704         \param name - Input - name of tag.
705         \param value - Input - Value to associate with given key.
706      */
707      ESCRIPT_DLL_API
708      void
709      setTaggedValueByName(std::string name,
710                           const boost::python::object& value);
711    
712      /**
713         \brief
714         Assign the given value to the tag. Implicitly converts this
715         object to type DataTagged if it is constant.
716    
717       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
718       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
719        ==>*
720    */    */
721      ESCRIPT_DLL_API
722    void    void
723    setTaggedValue(int tagKey,    setTaggedValue(int tagKey,
724                   const boost::python::object& value);                   const boost::python::object& value);
# Line 455  class Data { Line 726  class Data {
726    /**    /**
727       \brief       \brief
728       Assign the given value to the tag. Implicitly converts this       Assign the given value to the tag. Implicitly converts this
729       object to type DataTagged. Throws an exception if this object       object to type DataTagged if it is constant.
730       cannot be converted to a DataTagged object.  
731       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
732         \param pointshape - Input - The shape of the value parameter
733       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
734       Note: removed for now - this version not needed, and breaks escript.cpp       \param dataOffset - Input - Offset of the begining of the point within the value parameter
735    */    */
736    /*    ESCRIPT_DLL_API
737    void    void
738    setTaggedValue(int tagKey,    setTaggedValueFromCPP(int tagKey,
739                   const DataArrayView& value);              const DataTypes::ShapeType& pointshape,
740    */                          const DataTypes::ValueType& value,
741                int dataOffset=0);
742    
743    
744    
745    /**    /**
746      \brief      \brief
747      Copy other Data object into this Data object where mask is positive.      Copy other Data object into this Data object where mask is positive.
748    */    */
749      ESCRIPT_DLL_API
750    void    void
751    copyWithMask(const Data& other,    copyWithMask(const Data& other,
752                 const Data& mask);                 const Data& mask);
# Line 481  class Data { Line 757  class Data {
757    
758    /**    /**
759       \brief       \brief
760         set all values to zero
761         *
762      */
763      ESCRIPT_DLL_API
764      void
765      setToZero();
766    
767      /**
768         \brief
769       Interpolates this onto the given functionspace and returns       Interpolates this onto the given functionspace and returns
770       the result as a Data object.       the result as a Data object.
771         *
772    */    */
773      ESCRIPT_DLL_API
774    Data    Data
775    interpolate(const FunctionSpace& functionspace) const;    interpolate(const FunctionSpace& functionspace) const;
776    
777    
778      ESCRIPT_DLL_API
779      Data
780      interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
781                           double undef, Data& B, double Bmin, double Bstep);
782    
783    
784      // This signature needs to be tuned
785      ESCRIPT_DLL_API
786      Data
787      interpolateFromTable(boost::python::object table, double Amin, double Astep,
788                           double undef, Data& B, double Bmin, double Bstep);
789    
790    /**    /**
791       \brief       \brief
792       Calculates the gradient of the data at the data points of functionspace.       Calculates the gradient of the data at the data points of functionspace.
793       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.
794         *
795    */    */
796      ESCRIPT_DLL_API
797    Data    Data
798    gradOn(const FunctionSpace& functionspace) const;    gradOn(const FunctionSpace& functionspace) const;
799    
800      ESCRIPT_DLL_API
801    Data    Data
802    grad() const;    grad() const;
803    
804    /**    /**
805       \brief      \brief
806       Calculate the integral over the function space domain.       Calculate the integral over the function space domain as a python tuple.
807      */
808      ESCRIPT_DLL_API
809      boost::python::object
810      integrateToTuple_const() const;
811    
812    
813      /**
814        \brief
815         Calculate the integral over the function space domain as a python tuple.
816    */    */
817    boost::python::numeric::array    ESCRIPT_DLL_API
818    integrate() const;    boost::python::object
819      integrateToTuple();
820    
821    
822    
823    /**    /**
824       \brief       \brief
825         Returns 1./ Data object
826         *
827      */
828      ESCRIPT_DLL_API
829      Data
830      oneOver() const;
831      /**
832         \brief
833       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.
834         *
835    */    */
836      ESCRIPT_DLL_API
837    Data    Data
838    wherePositive() const;    wherePositive() const;
839    
840    /**    /**
841       \brief       \brief
842       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.
843         *
844    */    */
845      ESCRIPT_DLL_API
846    Data    Data
847    whereNegative() const;    whereNegative() const;
848    
849    /**    /**
850       \brief       \brief
851       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.
852         *
853    */    */
854      ESCRIPT_DLL_API
855    Data    Data
856    whereNonNegative() const;    whereNonNegative() const;
857    
858    /**    /**
859       \brief       \brief
860       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.
861         *
862    */    */
863      ESCRIPT_DLL_API
864    Data    Data
865    whereNonPositive() const;    whereNonPositive() const;
866    
867    /**    /**
868       \brief       \brief
869       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.
870         *
871    */    */
872      ESCRIPT_DLL_API
873    Data    Data
874    whereZero() const;    whereZero(double tol=0.0) const;
875    
876    /**    /**
877       \brief       \brief
878       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.
879         *
880      */
881      ESCRIPT_DLL_API
882      Data
883      whereNonZero(double tol=0.0) const;
884    
885      /**
886         \brief
887         Return the maximum absolute value of this Data object.
888    
889         The method is not const because lazy data needs to be expanded before Lsup can be computed.
890         The _const form can be used when the Data object is const, however this will only work for
891         Data which is not Lazy.
892    
893         For Data which contain no samples (or tagged Data for which no tags in use have a value)
894         zero is returned.
895      */
896      ESCRIPT_DLL_API
897      double
898      Lsup();
899    
900      ESCRIPT_DLL_API
901      double
902      Lsup_const() const;
903    
904    
905      /**
906         \brief
907         Return the maximum value of this Data object.
908    
909         The method is not const because lazy data needs to be expanded before sup can be computed.
910         The _const form can be used when the Data object is const, however this will only work for
911         Data which is not Lazy.
912    
913         For Data which contain no samples (or tagged Data for which no tags in use have a value)
914         a large negative value is returned.
915      */
916      ESCRIPT_DLL_API
917      double
918      sup();
919    
920      ESCRIPT_DLL_API
921      double
922      sup_const() const;
923    
924    
925      /**
926         \brief
927         Return the minimum value of this Data object.
928    
929         The method is not const because lazy data needs to be expanded before inf can be computed.
930         The _const form can be used when the Data object is const, however this will only work for
931         Data which is not Lazy.
932    
933         For Data which contain no samples (or tagged Data for which no tags in use have a value)
934         a large positive value is returned.
935      */
936      ESCRIPT_DLL_API
937      double
938      inf();
939    
940      ESCRIPT_DLL_API
941      double
942      inf_const() const;
943    
944    
945    
946      /**
947         \brief
948         Return the absolute value of each data point of this Data object.
949         *
950      */
951      ESCRIPT_DLL_API
952      Data
953      abs() const;
954    
955      /**
956         \brief
957         Return the maximum value of each data point of this Data object.
958         *
959      */
960      ESCRIPT_DLL_API
961      Data
962      maxval() const;
963    
964      /**
965         \brief
966         Return the minimum value of each data point of this Data object.
967         *
968      */
969      ESCRIPT_DLL_API
970      Data
971      minval() const;
972    
973      /**
974         \brief
975         Return the (sample number, data-point number) of the data point with
976         the minimum component value in this Data object.
977         \note If you are working in python, please consider using Locator
978    instead of manually manipulating process and point IDs.
979      */
980      ESCRIPT_DLL_API
981      const boost::python::tuple
982      minGlobalDataPoint() const;
983    
984      /**
985         \brief
986         Return the (sample number, data-point number) of the data point with
987         the minimum component value in this Data object.
988         \note If you are working in python, please consider using Locator
989    instead of manually manipulating process and point IDs.
990      */
991      ESCRIPT_DLL_API
992      const boost::python::tuple
993      maxGlobalDataPoint() const;
994    
995    
996    
997      /**
998         \brief
999         Return the sign of each data point of this Data object.
1000         -1 for negative values, zero for zero values, 1 for positive values.
1001         *
1002      */
1003      ESCRIPT_DLL_API
1004      Data
1005      sign() const;
1006    
1007      /**
1008         \brief
1009         Return the symmetric part of a matrix which is half the matrix plus its transpose.
1010         *
1011      */
1012      ESCRIPT_DLL_API
1013      Data
1014      symmetric() const;
1015    
1016      /**
1017         \brief
1018         Return the nonsymmetric part of a matrix which is half the matrix minus its transpose.
1019         *
1020      */
1021      ESCRIPT_DLL_API
1022      Data
1023      nonsymmetric() const;
1024    
1025      /**
1026         \brief
1027         Return the trace of a matrix
1028         *
1029      */
1030      ESCRIPT_DLL_API
1031      Data
1032      trace(int axis_offset) const;
1033    
1034      /**
1035         \brief
1036         Transpose each data point of this Data object around the given axis.
1037         *
1038      */
1039      ESCRIPT_DLL_API
1040      Data
1041      transpose(int axis_offset) const;
1042    
1043      /**
1044         \brief
1045         Return the eigenvalues of the symmetric part at each data point of this Data object in increasing values.
1046         Currently this function is restricted to rank 2, square shape, and dimension 3.
1047         *
1048      */
1049      ESCRIPT_DLL_API
1050      Data
1051      eigenvalues() const;
1052    
1053      /**
1054         \brief
1055         Return the eigenvalues and corresponding eigenvcetors of the symmetric part at each data point of this Data object.
1056         the eigenvalues are ordered in increasing size where eigenvalues with relative difference less than
1057         tol are treated as equal. The eigenvectors are orthogonal, normalized and the sclaed such that the
1058         first non-zero entry is positive.
1059         Currently this function is restricted to rank 2, square shape, and dimension 3
1060         *
1061      */
1062      ESCRIPT_DLL_API
1063      const boost::python::tuple
1064      eigenvalues_and_eigenvectors(const double tol=1.e-12) const;
1065    
1066      /**
1067         \brief
1068         swaps the components axis0 and axis1
1069         *
1070      */
1071      ESCRIPT_DLL_API
1072      Data
1073      swapaxes(const int axis0, const int axis1) const;
1074    
1075      /**
1076         \brief
1077         Return the error function erf of each data point of this Data object.
1078         *
1079    */    */
1080      ESCRIPT_DLL_API
1081    Data    Data
1082    whereNonZero() const;    erf() const;
1083    
1084    /**    /**
1085       \brief       \brief
1086       Return the sin of each data point of this Data object.       Return the sin of each data point of this Data object.
1087         *
1088    */    */
1089      ESCRIPT_DLL_API
1090    Data    Data
1091    sin() const;    sin() const;
1092    
1093    /**    /**
1094       \brief       \brief
1095       Return the cos of each data point of this Data object.       Return the cos of each data point of this Data object.
1096         *
1097    */    */
1098      ESCRIPT_DLL_API
1099    Data    Data
1100    cos() const;    cos() const;
1101    
1102    /**    /**
1103       \brief       \brief
1104       Return the tan of each data point of this Data object.       Return the tan of each data point of this Data object.
1105         *
1106    */    */
1107      ESCRIPT_DLL_API
1108    Data    Data
1109    tan() const;    tan() const;
1110    
1111    /**    /**
1112       \brief       \brief
1113       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.
1114         *
1115    */    */
1116      ESCRIPT_DLL_API
1117    Data    Data
1118    log() const;    asin() const;
1119    
1120    /**    /**
1121       \brief       \brief
1122       Return the natural log of each data point of this Data object.       Return the acos of each data point of this Data object.
1123         *
1124    */    */
1125      ESCRIPT_DLL_API
1126    Data    Data
1127    ln() const;    acos() const;
1128    
1129    /**    /**
1130       \brief       \brief
1131       Return the maximum absolute value of this Data object.       Return the atan of each data point of this Data object.
1132         *
1133    */    */
1134    double    ESCRIPT_DLL_API
1135    Lsup() const;    Data
1136      atan() const;
1137    
1138    /**    /**
1139       \brief       \brief
1140       Return the maximum value of this Data object.       Return the sinh of each data point of this Data object.
1141         *
1142    */    */
1143    double    ESCRIPT_DLL_API
1144    sup() const;    Data
1145      sinh() const;
1146    
1147    /**    /**
1148       \brief       \brief
1149       Return the minimum value of this Data object.       Return the cosh of each data point of this Data object.
1150         *
1151    */    */
1152    double    ESCRIPT_DLL_API
1153    inf() const;    Data
1154      cosh() const;
1155    
1156    /**    /**
1157       \brief       \brief
1158       Return the absolute value of each data point of this Data object.       Return the tanh of each data point of this Data object.
1159         *
1160    */    */
1161      ESCRIPT_DLL_API
1162    Data    Data
1163    abs() const;    tanh() const;
1164    
1165    /**    /**
1166       \brief       \brief
1167       Return the maximum value of each data point of this Data object.       Return the asinh of each data point of this Data object.
1168         *
1169    */    */
1170      ESCRIPT_DLL_API
1171    Data    Data
1172    maxval() const;    asinh() const;
1173    
1174    /**    /**
1175       \brief       \brief
1176       Return the minimum value of each data point of this Data object.       Return the acosh of each data point of this Data object.
1177         *
1178    */    */
1179      ESCRIPT_DLL_API
1180    Data    Data
1181    minval() const;    acosh() const;
1182    
1183    /**    /**
1184       \brief       \brief
1185       Return the length of each data point of this Data object.       Return the atanh of each data point of this Data object.
1186       sqrt(sum(A[i,j,k,l]^2))       *
1187    */    */
1188      ESCRIPT_DLL_API
1189    Data    Data
1190    length() const;    atanh() const;
1191    
1192    /**    /**
1193       \brief       \brief
1194       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.
1195       -1 for negative values, zero for zero values, 1 for positive values.       *
1196    */    */
1197      ESCRIPT_DLL_API
1198    Data    Data
1199    sign() const;    log10() const;
1200    
1201    /**    /**
1202      \transpose       \brief
1203      Transpose each data point of this Data object around the given axis.       Return the natural log of each data point of this Data object.
1204         *
1205    */    */
1206      ESCRIPT_DLL_API
1207    Data    Data
1208    transpose(int axis) const;    log() const;
1209    
1210    /**    /**
1211      \trace       \brief
1212      Calculate the trace of each data point of this Data object.       Return the exponential function of each data point of this Data object.
1213      sum(A[i,i,i,i])       *
1214    */    */
1215      ESCRIPT_DLL_API
1216    Data    Data
1217    trace() const;    exp() const;
1218    
1219    /**    /**
1220      \exp       \brief
1221      Return the exponential function of each data point of this Data object.       Return the square root of each data point of this Data object.
1222         *
1223    */    */
1224      ESCRIPT_DLL_API
1225    Data    Data
1226    exp() const;    sqrt() const;
1227    
1228    /**    /**
1229      \sqrt       \brief
1230      Return the square root of each data point of this Data object.       Return the negation of each data point of this Data object.
1231         *
1232    */    */
1233      ESCRIPT_DLL_API
1234    Data    Data
1235    sqrt() const;    neg() const;
1236    
1237      /**
1238         \brief
1239         Return the identity of each data point of this Data object.
1240         Simply returns this object unmodified.
1241         *
1242      */
1243      ESCRIPT_DLL_API
1244      Data
1245      pos() const;
1246    
1247    /**    /**
1248       \brief       \brief
1249       Return the given power of each data point of this Data object.       Return the given power of each data point of this Data object.
1250    
1251         \param right Input - the power to raise the object to.
1252         *
1253    */    */
1254      ESCRIPT_DLL_API
1255    Data    Data
1256    powD(const Data& right) const;    powD(const Data& right) const;
1257    
1258      /**
1259         \brief
1260         Return the given power of each data point of this boost python object.
1261    
1262         \param right Input - the power to raise the object to.
1263         *
1264       */
1265      ESCRIPT_DLL_API
1266    Data    Data
1267    powO(const boost::python::object& right) const;    powO(const boost::python::object& right) const;
1268    
1269    /**    /**
1270      \brief       \brief
1271      Return the negation of each data point of this Data object.       Return the given power of each data point of this boost python object.
1272    */  
1273         \param left Input - the bases
1274         *
1275       */
1276    
1277      ESCRIPT_DLL_API
1278    Data    Data
1279    neg() const;    rpowO(const boost::python::object& left) const;
1280    
1281    /**    /**
1282      \brief       \brief
1283      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.  
1284    */    */
1285    Data    ESCRIPT_DLL_API
1286    pos() const;    void
1287      saveDX(std::string fileName) const;
1288    
1289      /**
1290         \brief
1291         writes the object to a file in the VTK file format
1292      */
1293      ESCRIPT_DLL_API
1294      void
1295      saveVTK(std::string fileName) const;
1296    
1297    
1298    
1299    /**    /**
1300       \brief       \brief
1301       Overloaded operator +=       Overloaded operator +=
1302       \param right - Input - The right hand side.       \param right - Input - The right hand side.
1303         *
1304    */    */
1305      ESCRIPT_DLL_API
1306    Data& operator+=(const Data& right);    Data& operator+=(const Data& right);
1307      ESCRIPT_DLL_API
1308    Data& operator+=(const boost::python::object& right);    Data& operator+=(const boost::python::object& right);
1309    
1310      ESCRIPT_DLL_API
1311      Data& operator=(const Data& other);
1312    
1313    /**    /**
1314       \brief       \brief
1315       Overloaded operator -=       Overloaded operator -=
1316       \param right - Input - The right hand side.       \param right - Input - The right hand side.
1317         *
1318    */    */
1319      ESCRIPT_DLL_API
1320    Data& operator-=(const Data& right);    Data& operator-=(const Data& right);
1321      ESCRIPT_DLL_API
1322    Data& operator-=(const boost::python::object& right);    Data& operator-=(const boost::python::object& right);
1323    
1324   /**   /**
1325       \brief       \brief
1326       Overloaded operator *=       Overloaded operator *=
1327       \param right - Input - The right hand side.       \param right - Input - The right hand side.
1328         *
1329    */    */
1330      ESCRIPT_DLL_API
1331    Data& operator*=(const Data& right);    Data& operator*=(const Data& right);
1332      ESCRIPT_DLL_API
1333    Data& operator*=(const boost::python::object& right);    Data& operator*=(const boost::python::object& right);
1334    
1335   /**   /**
1336       \brief       \brief
1337       Overloaded operator /=       Overloaded operator /=
1338       \param right - Input - The right hand side.       \param right - Input - The right hand side.
1339         *
1340    */    */
1341      ESCRIPT_DLL_API
1342    Data& operator/=(const Data& right);    Data& operator/=(const Data& right);
1343      ESCRIPT_DLL_API
1344    Data& operator/=(const boost::python::object& right);    Data& operator/=(const boost::python::object& right);
1345    
1346    /**    /**
1347       \brief       \brief
1348       Returns true if this can be interpolated to functionspace.       Returns true if this can be interpolated to functionspace.
1349    */    */
1350      ESCRIPT_DLL_API
1351    bool    bool
1352    probeInterpolation(const FunctionSpace& functionspace) const;    probeInterpolation(const FunctionSpace& functionspace) const;
1353    
# Line 748  class Data { Line 1366  class Data {
1366       \param key - Input - python slice tuple specifying       \param key - Input - python slice tuple specifying
1367       slice to return.       slice to return.
1368    */    */
1369      ESCRIPT_DLL_API
1370    Data    Data
1371    getItem(const boost::python::object& key) const;    getItem(const boost::python::object& key) const;
1372    
# Line 755  class Data { Line 1374  class Data {
1374       \brief       \brief
1375       Copies slice from value into this Data object.       Copies slice from value into this Data object.
1376    
      \description  
1377       Implements the [] set operator in python.       Implements the [] set operator in python.
1378       Calls setSlice.       Calls setSlice.
1379    
# Line 763  class Data { Line 1381  class Data {
1381       slice to copy from value.       slice to copy from value.
1382       \param value - Input - Data object to copy from.       \param value - Input - Data object to copy from.
1383    */    */
1384      ESCRIPT_DLL_API
1385    void    void
1386    setItemD(const boost::python::object& key,    setItemD(const boost::python::object& key,
1387             const Data& value);             const Data& value);
1388    
1389      ESCRIPT_DLL_API
1390    void    void
1391    setItemO(const boost::python::object& key,    setItemO(const boost::python::object& key,
1392             const boost::python::object& value);             const boost::python::object& value);
# Line 779  class Data { Line 1399  class Data {
1399       this Data object.       this Data object.
1400    */    */
1401    template <class UnaryFunction>    template <class UnaryFunction>
1402      ESCRIPT_DLL_API
1403    inline    inline
1404    void    void
1405    unaryOp(UnaryFunction operation);    unaryOp2(UnaryFunction operation);
1406    
1407    /**    /**
1408       \brief       \brief
1409       Return a Data object containing the specified slice of       Return a Data object containing the specified slice of
1410       this Data object.       this Data object.
1411       \param region - Input - Region to copy.       \param region - Input - Region to copy.
1412         *
1413    */    */
1414      ESCRIPT_DLL_API
1415    Data    Data
1416    getSlice(const DataArrayView::RegionType& region) const;    getSlice(const DataTypes::RegionType& region) const;
1417    
1418    /**    /**
1419       \brief       \brief
# Line 798  class Data { Line 1421  class Data {
1421       Data object.       Data object.
1422       \param value - Input - Data to copy from.       \param value - Input - Data to copy from.
1423       \param region - Input - Region to copy.       \param region - Input - Region to copy.
1424         *
1425    */    */
1426      ESCRIPT_DLL_API
1427    void    void
1428    setSlice(const Data& value,    setSlice(const Data& value,
1429             const DataArrayView::RegionType& region);             const DataTypes::RegionType& region);
1430    
1431      /**
1432         \brief
1433         print the data values to stdout. Used for debugging
1434      */
1435      ESCRIPT_DLL_API
1436      void
1437            print(void);
1438    
1439      /**
1440         \brief
1441         return the MPI rank number of the local data
1442                     MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()
1443                     is returned
1444      */
1445      ESCRIPT_DLL_API
1446            int
1447            get_MPIRank(void) const;
1448    
1449      /**
1450         \brief
1451         return the MPI rank number of the local data
1452                     MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()
1453                     is returned
1454      */
1455      ESCRIPT_DLL_API
1456            int
1457            get_MPISize(void) const;
1458    
1459      /**
1460         \brief
1461         return the MPI rank number of the local data
1462                     MPI_COMM_WORLD is assumed and returned.
1463      */
1464      ESCRIPT_DLL_API
1465            MPI_Comm
1466            get_MPIComm(void) const;
1467    
1468      /**
1469         \brief
1470         return the object produced by the factory, which is a DataConstant or DataExpanded
1471        TODO Ownership of this object should be explained in doco.
1472      */
1473      ESCRIPT_DLL_API
1474            DataAbstract*
1475            borrowData(void) const;
1476    
1477      ESCRIPT_DLL_API
1478            DataAbstract_ptr
1479            borrowDataPtr(void) const;
1480    
1481      ESCRIPT_DLL_API
1482            DataReady_ptr
1483            borrowReadyPtr(void) const;
1484    
1485    
1486    
1487      /**
1488         \brief
1489         Return a pointer to the beginning of the datapoint at the specified offset.
1490         TODO Eventually these should be inlined.
1491         \param i - position(offset) in the underlying datastructure
1492      */
1493    
1494      ESCRIPT_DLL_API
1495            DataTypes::ValueType::const_reference
1496            getDataAtOffsetRO(DataTypes::ValueType::size_type i);
1497    
1498    
1499      ESCRIPT_DLL_API
1500            DataTypes::ValueType::reference
1501            getDataAtOffsetRW(DataTypes::ValueType::size_type i);
1502    
1503    
1504    
1505    /**
1506       \brief Create a buffer for use by getSample
1507       Allocates a DataVector large enough for DataLazy::resolveSample to operate on for the current Data.
1508       Do not use this buffer for other Data instances (unless you are sure they will be the same size).
1509      
1510       In multi-threaded sections, this needs to be called on each thread.
1511    
1512       \return A BufferGroup* if Data is lazy, NULL otherwise.
1513       \warning This pointer must be deallocated using freeSampleBuffer to avoid cross library memory issues.
1514    */
1515      ESCRIPT_DLL_API
1516      BufferGroup*
1517      allocSampleBuffer() const;
1518    
1519    /**
1520       \brief Free a buffer allocated with allocSampleBuffer.
1521       \param buffer Input - pointer to the buffer to deallocate.
1522    */
1523    ESCRIPT_DLL_API void freeSampleBuffer(BufferGroup* buffer);
1524    
1525   protected:   protected:
1526    
1527   private:   private:
1528    
1529      double
1530      LsupWorker() const;
1531    
1532      double
1533      supWorker() const;
1534    
1535      double
1536      infWorker() const;
1537    
1538      boost::python::object
1539      integrateWorker() const;
1540    
1541      void
1542      calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1543    
1544      void
1545      calc_maxGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
1546    
1547    
1548    /**    /**
1549       \brief       \brief
1550       Check *this and the right operand are compatible. Throws       Check *this and the right operand are compatible. Throws
# Line 823  class Data { Line 1561  class Data {
1561    /**    /**
1562       \brief       \brief
1563       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
1564       this Data object and return the single double value result.       this Data object according to the given function and return the single value result.
1565    */    */
1566    template <class UnaryFunction>    template <class BinaryFunction>
1567    inline    inline
1568    double    double
1569    algorithm(UnaryFunction operation) const;    algorithm(BinaryFunction operation,
1570                double initial_value) const;
1571    
1572    /**    /**
1573       \brief       \brief
1574       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
1575       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
1576       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
1577       is expanded *this will be expanded if necessary.       this Data object
      RHS is a Data object.  
1578    */    */
1579    template <class BinaryFunction>    template <class BinaryFunction>
1580    inline    inline
1581    void    Data
1582    binaryOp(const Data& right,    dp_algorithm(BinaryFunction operation,
1583             BinaryFunction operation);                 double initial_value) const;
1584    
1585    /**    /**
1586       \brief       \brief
1587       Perform the given binary operation on all of the data's elements.       Perform the given binary operation on all of the data's elements.
1588       RHS is a boost::python object.       The underlying type of the right hand side (right) determines the final
1589         type of *this after the operation. For example if the right hand side
1590         is expanded *this will be expanded if necessary.
1591         RHS is a Data object.
1592    */    */
1593    template <class BinaryFunction>    template <class BinaryFunction>
1594    inline    inline
1595    void    void
1596    binaryOp(const boost::python::object& right,    binaryOp(const Data& right,
1597             BinaryFunction operation);             BinaryFunction operation);
1598    
1599    /**    /**
# Line 875  class Data { Line 1616  class Data {
1616       \brief       \brief
1617       Construct a Data object of the appropriate type.       Construct a Data object of the appropriate type.
1618    */    */
1619    template <class IValueType>  
1620    void    void
1621    initialise(const IValueType& value,    initialise(const DataTypes::ValueType& value,
1622             const DataTypes::ShapeType& shape,
1623               const FunctionSpace& what,               const FunctionSpace& what,
1624               bool expanded);               bool expanded);
1625    
   /**  
      \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.  
   */  
1626    void    void
1627    reshapeDataPoint(const DataArrayView::ShapeType& shape);    initialise(const WrappedArray& value,
1628                     const FunctionSpace& what,
1629                     bool expanded);
1630    
1631      //
1632      // flag to protect the data object against any update
1633      bool m_protected;
1634      mutable bool m_shared;
1635      bool m_lazy;
1636    
1637    //    //
1638    // pointer to the actual data object    // pointer to the actual data object
1639    boost::shared_ptr<DataAbstract> m_data;  //   boost::shared_ptr<DataAbstract> m_data;
1640      DataAbstract_ptr m_data;
1641    
1642  };  // If possible please use getReadyPtr instead.
1643    // But see warning below.
1644      const DataReady*
1645      getReady() const;
1646    
1647  template <class IValueType>    DataReady*
1648  void    getReady();
1649  Data::initialise(const IValueType& value,  
1650                   const FunctionSpace& what,  
1651                   bool expanded)  // Be wary of using this for local operations since it (temporarily) increases reference count.
1652  {  // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1653    // getReady() instead
1654      DataReady_ptr
1655      getReadyPtr();
1656    
1657      const_DataReady_ptr
1658      getReadyPtr() const;
1659    
1660    
1661      /**
1662       \brief Update the Data's shared flag
1663       This indicates that the DataAbstract used by this object is now shared (or no longer shared).
1664       For internal use only.
1665      */
1666      void updateShareStatus(bool nowshared) const
1667      {
1668        m_shared=nowshared;     // m_shared is mutable
1669      }
1670    
1671      // In the isShared() method below:
1672      // A problem would occur if m_data (the address pointed to) were being modified
1673      // while the call m_data->is_shared is being executed.
1674      //
1675      // Q: So why do I think this code can be thread safe/correct?
1676      // A: We need to make some assumptions.
1677      // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1678      // 2. We assume that no constructions or assignments which will share previously unshared
1679      //    will occur while this call is executing. This is consistent with the way Data:: and C are written.
1680    //    //
1681    // 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
1682    // 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.
1683    // 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.
1684    // within the shared_ptr constructor.    // For any threads executing before the flag switches they will assume the object is still shared.
1685    if (expanded) {    bool isShared() const
1686      DataAbstract* temp=new DataExpanded(value,what);    {
1687      boost::shared_ptr<DataAbstract> temp_data(temp);      return m_shared;
1688      m_data=temp_data;  /*  if (m_shared) return true;
1689    } else {      if (m_data->isShared())        
1690      DataAbstract* temp=new DataConstant(value,what);      {                  
1691      boost::shared_ptr<DataAbstract> temp_data(temp);          updateShareStatus(true);
1692      m_data=temp_data;          return true;
1693        }
1694        return false;*/
1695      }
1696    
1697      void forceResolve()
1698      {
1699        if (isLazy())
1700        {
1701            #ifdef _OPENMP
1702            if (omp_in_parallel())
1703            {   // Yes this is throwing an exception out of an omp thread which is forbidden.
1704            throw DataException("Please do not call forceResolve() in a parallel region.");
1705            }
1706            #endif
1707            resolve();
1708        }
1709      }
1710    
1711      /**
1712      \brief if another object is sharing out member data make a copy to work with instead.
1713      This code should only be called from single threaded sections of code.
1714      */
1715      void exclusiveWrite()
1716      {
1717    #ifdef _OPENMP
1718      if (omp_in_parallel())
1719      {
1720    // *((int*)0)=17;
1721        throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1722      }
1723    #endif
1724        forceResolve();
1725        if (isShared())
1726        {
1727            DataAbstract* t=m_data->deepCopy();
1728            set_m_data(DataAbstract_ptr(t));
1729        }
1730      }
1731    
1732      /**
1733      \brief checks if caller can have exclusive write to the object
1734      */
1735      void checkExclusiveWrite()
1736      {
1737        if  (isLazy() || isShared())
1738        {
1739            throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
1740        }
1741    }    }
1742    
1743      /**
1744      \brief Modify the data abstract hosted by this Data object
1745      For internal use only.
1746      Passing a pointer to null is permitted (do this in the destructor)
1747      \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1748      */
1749      void set_m_data(DataAbstract_ptr p);
1750    
1751      friend class DataAbstract;        // To allow calls to updateShareStatus
1752    
1753    };
1754    
1755    }   // end namespace escript
1756    
1757    
1758    // No, this is not supposed to be at the top of the file
1759    // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1760    // so that I can dynamic cast between them below.
1761    #include "DataReady.h"
1762    #include "DataLazy.h"
1763    
1764    namespace escript
1765    {
1766    
1767    inline
1768    const DataReady*
1769    Data::getReady() const
1770    {
1771       const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1772       EsysAssert((dr!=0), "Error - casting to DataReady.");
1773       return dr;
1774  }  }
1775    
1776    inline
1777    DataReady*
1778    Data::getReady()
1779    {
1780       DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1781       EsysAssert((dr!=0), "Error - casting to DataReady.");
1782       return dr;
1783    }
1784    
1785    // Be wary of using this for local operations since it (temporarily) increases reference count.
1786    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1787    // getReady() instead
1788    inline
1789    DataReady_ptr
1790    Data::getReadyPtr()
1791    {
1792       DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1793       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1794       return dr;
1795    }
1796    
1797    
1798    inline
1799    const_DataReady_ptr
1800    Data::getReadyPtr() const
1801    {
1802       const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1803       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1804       return dr;
1805    }
1806    
1807    inline
1808    DataAbstract::ValueType::value_type*
1809    Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1810    {
1811       if (isLazy())
1812       {
1813        throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1814       }
1815       return getReady()->getSampleDataRW(sampleNo);
1816    }
1817    
1818    inline
1819    const DataAbstract::ValueType::value_type*
1820    Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg)
1821    {
1822       DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1823       if (l!=0)
1824       {
1825        size_t offset=0;
1826        if (bufferg==NULL)
1827        {
1828            throw DataException("Error, attempt to getSampleDataRO for lazy Data with buffer==NULL");
1829        }
1830        const DataTypes::ValueType* res=l->resolveSample(*bufferg,sampleNo,offset);
1831        return &((*res)[offset]);
1832       }
1833       return getReady()->getSampleDataRO(sampleNo);
1834    }
1835    
1836    
1837    
1838    /**
1839       Modify a filename for MPI parallel output to multiple files
1840    */
1841    char *Escript_MPI_appendRankToFileName(const char *, int, int);
1842    
1843  /**  /**
1844     Binary Data object operators.     Binary Data object operators.
1845  */  */
1846    inline double rpow(double x,double y)
1847    {
1848        return pow(y,x);
1849    }
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
1860    Operator-    Operator-
1861    Takes two Data objects.    Takes two Data objects.
1862  */  */
1863  Data operator-(const Data& left, const Data& right);  ESCRIPT_DLL_API Data operator-(const Data& left, const Data& right);
1864    
1865  /**  /**
1866    \brief    \brief
1867    Operator*    Operator*
1868    Takes two Data objects.    Takes two Data objects.
1869  */  */
1870  Data operator*(const Data& left, const Data& right);  ESCRIPT_DLL_API Data operator*(const Data& left, const Data& right);
1871    
1872  /**  /**
1873    \brief    \brief
1874    Operator/    Operator/
1875    Takes two Data objects.    Takes two Data objects.
1876  */  */
1877  Data operator/(const Data& left, const Data& right);  ESCRIPT_DLL_API Data operator/(const Data& left, const Data& right);
1878    
1879  /**  /**
1880    \brief    \brief
# Line 957  Data operator/(const Data& left, const D Line 1882  Data operator/(const Data& left, const D
1882    Takes LHS Data object and RHS python::object.    Takes LHS Data object and RHS python::object.
1883    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1884  */  */
1885  Data operator+(const Data& left, const boost::python::object& right);  ESCRIPT_DLL_API Data operator+(const Data& left, const boost::python::object& right);
1886    
1887  /**  /**
1888    \brief    \brief
# Line 965  Data operator+(const Data& left, const b Line 1890  Data operator+(const Data& left, const b
1890    Takes LHS Data object and RHS python::object.    Takes LHS Data object and RHS python::object.
1891    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1892  */  */
1893  Data operator-(const Data& left, const boost::python::object& right);  ESCRIPT_DLL_API Data operator-(const Data& left, const boost::python::object& right);
1894    
1895  /**  /**
1896    \brief    \brief
# Line 973  Data operator-(const Data& left, const b Line 1898  Data operator-(const Data& left, const b
1898    Takes LHS Data object and RHS python::object.    Takes LHS Data object and RHS python::object.
1899    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1900  */  */
1901  Data operator*(const Data& left, const boost::python::object& right);  ESCRIPT_DLL_API Data operator*(const Data& left, const boost::python::object& right);
1902    
1903  /**  /**
1904    \brief    \brief
# Line 981  Data operator*(const Data& left, const b Line 1906  Data operator*(const Data& left, const b
1906    Takes LHS Data object and RHS python::object.    Takes LHS Data object and RHS python::object.
1907    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1908  */  */
1909  Data operator/(const Data& left, const boost::python::object& right);  ESCRIPT_DLL_API Data operator/(const Data& left, const boost::python::object& right);
1910    
1911  /**  /**
1912    \brief    \brief
# Line 989  Data operator/(const Data& left, const b Line 1914  Data operator/(const Data& left, const b
1914    Takes LHS python::object and RHS Data object.    Takes LHS python::object and RHS Data object.
1915    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1916  */  */
1917  Data operator+(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator+(const boost::python::object& left, const Data& right);
1918    
1919  /**  /**
1920    \brief    \brief
# Line 997  Data operator+(const boost::python::obje Line 1922  Data operator+(const boost::python::obje
1922    Takes LHS python::object and RHS Data object.    Takes LHS python::object and RHS Data object.
1923    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1924  */  */
1925  Data operator-(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator-(const boost::python::object& left, const Data& right);
1926    
1927  /**  /**
1928    \brief    \brief
# Line 1005  Data operator-(const boost::python::obje Line 1930  Data operator-(const boost::python::obje
1930    Takes LHS python::object and RHS Data object.    Takes LHS python::object and RHS Data object.
1931    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1932  */  */
1933  Data operator*(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator*(const boost::python::object& left, const Data& right);
1934    
1935  /**  /**
1936    \brief    \brief
# Line 1013  Data operator*(const boost::python::obje Line 1938  Data operator*(const boost::python::obje
1938    Takes LHS python::object and RHS Data object.    Takes LHS python::object and RHS Data object.
1939    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1940  */  */
1941  Data operator/(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator/(const boost::python::object& left, const Data& right);
1942    
1943    
1944    
1945  /**  /**
1946    \brief    \brief
1947    Output operator    Output operator
1948  */  */
1949  std::ostream& operator<<(std::ostream& o, const Data& data);  ESCRIPT_DLL_API std::ostream& operator<<(std::ostream& o, const Data& data);
1950    
1951  /**  /**
1952    \brief    \brief
1953    Return true if operands are equivalent, else return false.    Compute a tensor product of two Data objects
1954    NB: this operator does very little at this point, and isn't to    \param arg_0 - Input - Data object
1955    be relied on. Requires further implementation.    \param arg_1 - Input - Data object
1956      \param axis_offset - Input - axis offset
1957      \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1958  */  */
1959  //bool operator==(const Data& left, const Data& right);  ESCRIPT_DLL_API
1960    Data
1961    C_GeneralTensorProduct(Data& arg_0,
1962                         Data& arg_1,
1963                         int axis_offset=0,
1964                         int transpose=0);
1965    
1966  /**  /**
1967    \brief    \brief
# Line 1042  Data::binaryOp(const Data& right, Line 1976  Data::binaryOp(const Data& right,
1976  {  {
1977     //     //
1978     // 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
1979     if (getPointDataView().getRank()==0 && right.getPointDataView().getRank()!=0) {     if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1980       reshapeDataPoint(right.getPointDataView().getShape());       throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
1981       }
1982    
1983       if (isLazy() || right.isLazy())
1984       {
1985         throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
1986     }     }
1987     //     //
1988     // initially make the temporary a shallow copy     // initially make the temporary a shallow copy
1989     Data tempRight(right);     Data tempRight(right);
1990    
1991     if (getFunctionSpace()!=right.getFunctionSpace()) {     if (getFunctionSpace()!=right.getFunctionSpace()) {
1992       if (right.probeInterpolation(getFunctionSpace())) {       if (right.probeInterpolation(getFunctionSpace())) {
1993         //         //
1994         // an interpolation is required so create a new Data         // an interpolation is required so create a new Data
1995         tempRight=Data(right,this->getFunctionSpace());         tempRight=Data(right,this->getFunctionSpace());
1996       } else if (probeInterpolation(right.getFunctionSpace())) {       } else if (probeInterpolation(right.getFunctionSpace())) {
1997         //         //
1998         // interpolate onto the RHS function space         // interpolate onto the RHS function space
1999         Data tempLeft(*this,right.getFunctionSpace());         Data tempLeft(*this,right.getFunctionSpace());
2000         m_data=tempLeft.m_data;  //        m_data=tempLeft.m_data;
2001           set_m_data(tempLeft.m_data);
2002       }       }
2003     }     }
2004     operandCheck(tempRight);     operandCheck(tempRight);
# Line 1073  Data::binaryOp(const Data& right, Line 2014  Data::binaryOp(const Data& right,
2014       // of any data type       // of any data type
2015       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2016       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2017       escript::binaryOp(*leftC,*(tempRight.m_data.get()),operation);       escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
2018     } else if (isTagged()) {     } else if (isTagged()) {
2019       //       //
2020       // 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 2040  Data::binaryOp(const Data& right,
2040    
2041  /**  /**
2042    \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  
2043    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.
2044    Given operation combines each element of each data point, thus argument    Given operation combines each element of each data point, thus argument
2045    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.
2046    Calls escript::algorithm.    Calls escript::algorithm.
2047  */  */
2048  template <class UnaryFunction>  template <class BinaryFunction>
2049  inline  inline
2050  double  double
2051  Data::algorithm(UnaryFunction operation) const  Data::algorithm(BinaryFunction operation, double initial_value) const
2052  {  {
2053    if (isExpanded()) {    if (isExpanded()) {
2054      DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2055      EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");      EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2056      return escript::algorithm(*leftC,operation);      return escript::algorithm(*leftC,operation,initial_value);
2057    } else if (isTagged()) {    } else if (isTagged()) {
2058      DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
2059      EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");      EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");
2060      return escript::algorithm(*leftC,operation);      return escript::algorithm(*leftC,operation,initial_value);
2061    } else if (isConstant()) {    } else if (isConstant()) {
2062      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2063      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2064      return escript::algorithm(*leftC,operation);      return escript::algorithm(*leftC,operation,initial_value);
2065      } else if (isEmpty()) {
2066        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2067      } else if (isLazy()) {
2068        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2069      } else {
2070        throw DataException("Error - Data encapsulates an unknown type.");
2071    }    }
   return 0;  
2072  }  }
2073    
2074  /**  /**
2075    \brief    \brief
2076    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.
2077    Given operation combines each element within each data point into a scalar,    Given operation combines each element within each data point into a scalar,
2078    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
2079    rank 0 Data object.    rank 0 Data object.
2080    Calls escript::dp_algorithm.    Calls escript::dp_algorithm.
2081  */  */
2082  template <class UnaryFunction>  template <class BinaryFunction>
2083  inline  inline
2084  Data  Data
2085  dp_algorithm(const Data& data,  Data::dp_algorithm(BinaryFunction operation, double initial_value) const
              UnaryFunction operation)  
2086  {  {
2087    Data result(0,DataArrayView::ShapeType(),data.getFunctionSpace(),data.isExpanded());    if (isEmpty()) {
2088    if (data.isExpanded()) {      throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2089      DataExpanded* dataE=dynamic_cast<DataExpanded*>(data.m_data.get());    }
2090      else if (isExpanded()) {
2091        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2092        DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2093      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2094      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2095      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2096      escript::dp_algorithm(*dataE,*resultE,operation);      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2097    } else if (data.isTagged()) {      return result;
2098      DataTagged* dataT=dynamic_cast<DataTagged*>(data.m_data.get());    }
2099      DataTagged* resultT=dynamic_cast<DataTagged*>(result.m_data.get());    else if (isTagged()) {
2100        DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());
2101      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2102      EsysAssert((resultT!=0), "Programming error - casting result to DataTagged.");      DataTypes::ValueType defval(1);
2103      escript::dp_algorithm(*dataT,*resultT,operation);      defval[0]=0;
2104    } else if (data.isConstant()) {      DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2105      DataConstant* dataC=dynamic_cast<DataConstant*>(data.m_data.get());      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2106        return Data(resultT);   // note: the Data object now owns the resultT pointer
2107      }
2108      else if (isConstant()) {
2109        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2110        DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2111      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2112      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2113      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2114      escript::dp_algorithm(*dataC,*resultC,operation);      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2115        return result;
2116      } else if (isLazy()) {
2117        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2118      } else {
2119        throw DataException("Error - Data encapsulates an unknown type.");
2120    }    }
2121    return result;  }
2122    
2123    /**
2124      \brief
2125      Compute a tensor operation with two Data objects
2126      \param arg_0 - Input - Data object
2127      \param arg_1 - Input - Data object
2128      \param operation - Input - Binary op functor
2129    */
2130    template <typename BinaryFunction>
2131    inline
2132    Data
2133    C_TensorBinaryOperation(Data const &arg_0,
2134                            Data const &arg_1,
2135                            BinaryFunction operation)
2136    {
2137      if (arg_0.isEmpty() || arg_1.isEmpty())
2138      {
2139         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2140      }
2141      if (arg_0.isLazy() || arg_1.isLazy())
2142      {
2143         throw DataException("Error - Operations not permitted on lazy data.");
2144      }
2145      // Interpolate if necessary and find an appropriate function space
2146      Data arg_0_Z, arg_1_Z;
2147      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2148        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2149          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2150          arg_1_Z = Data(arg_1);
2151        }
2152        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2153          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2154          arg_0_Z =Data(arg_0);
2155        }
2156        else {
2157          throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible function spaces.");
2158        }
2159      } else {
2160          arg_0_Z = Data(arg_0);
2161          arg_1_Z = Data(arg_1);
2162      }
2163      // Get rank and shape of inputs
2164      int rank0 = arg_0_Z.getDataPointRank();
2165      int rank1 = arg_1_Z.getDataPointRank();
2166      DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2167      DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2168      int size0 = arg_0_Z.getDataPointSize();
2169      int size1 = arg_1_Z.getDataPointSize();
2170      // Declare output Data object
2171      Data res;
2172    
2173      if (shape0 == shape1) {
2174        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2175          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2176          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2177          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2178          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2179    
2180          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2181        }
2182        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2183    
2184          // Prepare the DataConstant input
2185          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2186    
2187          // Borrow DataTagged input from Data object
2188          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2189    
2190          // Prepare a DataTagged output 2
2191          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2192          res.tag();
2193          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2194    
2195          // Prepare offset into DataConstant
2196          int offset_0 = tmp_0->getPointOffset(0,0);
2197          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2198    
2199          // Get the pointers to the actual data
2200          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2201          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2202    
2203          // Compute a result for the default
2204          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2205          // Compute a result for each tag
2206          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2207          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2208          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2209            tmp_2->addTag(i->first);
2210            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2211            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2212    
2213            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2214          }
2215    
2216        }
2217        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2218          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2219          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2220          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2221          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2222    
2223          int sampleNo_1,dataPointNo_1;
2224          int numSamples_1 = arg_1_Z.getNumSamples();
2225          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2226          int offset_0 = tmp_0->getPointOffset(0,0);
2227          res.requireWrite();
2228          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2229          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2230            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2231              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2232              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2233              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2234              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2235              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2236              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2237            }
2238          }
2239    
2240        }
2241        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2242          // Borrow DataTagged input from Data object
2243          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2244    
2245          // Prepare the DataConstant input
2246          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2247    
2248          // Prepare a DataTagged output 2
2249          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2250          res.tag();
2251          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2252    
2253          // Prepare offset into DataConstant
2254          int offset_1 = tmp_1->getPointOffset(0,0);
2255    
2256          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2257          // Get the pointers to the actual data
2258          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2259          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2260          // Compute a result for the default
2261          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2262          // Compute a result for each tag
2263          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2264          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2265          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2266            tmp_2->addTag(i->first);
2267            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2268            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2269            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2270          }
2271    
2272        }
2273        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2274          // Borrow DataTagged input from Data object
2275          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2276    
2277          // Borrow DataTagged input from Data object
2278          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2279    
2280          // Prepare a DataTagged output 2
2281          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2282          res.tag();        // DataTagged output
2283          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2284    
2285          // Get the pointers to the actual data
2286          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2287          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2288          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2289    
2290          // Compute a result for the default
2291          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2292          // Merge the tags
2293          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2294          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2295          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2296          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2297            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2298          }
2299          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2300            tmp_2->addTag(i->first);
2301          }
2302          // Compute a result for each tag
2303          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2304          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2305    
2306            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2307            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2308            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2309    
2310            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2311          }
2312    
2313        }
2314        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2315          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2316          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2317          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2318          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2319          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2320    
2321          int sampleNo_0,dataPointNo_0;
2322          int numSamples_0 = arg_0_Z.getNumSamples();
2323          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2324          res.requireWrite();
2325          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2326          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2327            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2328            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2329            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2330              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2331              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2332              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2333              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2334              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2335            }
2336          }
2337    
2338        }
2339        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2340          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2341          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2342          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2343          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2344    
2345          int sampleNo_0,dataPointNo_0;
2346          int numSamples_0 = arg_0_Z.getNumSamples();
2347          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2348          int offset_1 = tmp_1->getPointOffset(0,0);
2349          res.requireWrite();
2350          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2351          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2352            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2353              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2354              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2355    
2356              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2357              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2358              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2359    
2360    
2361              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2362            }
2363          }
2364    
2365        }
2366        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2367          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2368          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2369          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2370          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2371          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2372    
2373          int sampleNo_0,dataPointNo_0;
2374          int numSamples_0 = arg_0_Z.getNumSamples();
2375          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2376          res.requireWrite();
2377          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2378          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2379            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2380            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2381            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2382              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2383              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2384              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2385              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2386              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2387            }
2388          }
2389    
2390        }
2391        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2392          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2393          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2394          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2395          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2396          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2397    
2398          int sampleNo_0,dataPointNo_0;
2399          int numSamples_0 = arg_0_Z.getNumSamples();
2400          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2401          res.requireWrite();
2402          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2403          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2404            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2405              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2406              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2407              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2408              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2409              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2410              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2411              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2412            }
2413          }
2414    
2415        }
2416        else {
2417          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2418        }
2419    
2420      } else if (0 == rank0) {
2421        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2422          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
2423          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2424          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2425          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2426          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2427        }
2428        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2429    
2430          // Prepare the DataConstant input
2431          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2432    
2433          // Borrow DataTagged input from Data object
2434          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2435    
2436          // Prepare a DataTagged output 2
2437          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataTagged output
2438          res.tag();
2439          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2440    
2441          // Prepare offset into DataConstant
2442          int offset_0 = tmp_0->getPointOffset(0,0);
2443          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2444    
2445          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2446          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2447    
2448          // Compute a result for the default
2449          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2450          // Compute a result for each tag
2451          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2452          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2453          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2454            tmp_2->addTag(i->first);
2455            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2456            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2457            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2458          }
2459    
2460        }
2461        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2462    
2463          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2464          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2465          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2466          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2467    
2468          int sampleNo_1,dataPointNo_1;
2469          int numSamples_1 = arg_1_Z.getNumSamples();
2470          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2471          int offset_0 = tmp_0->getPointOffset(0,0);
2472          res.requireWrite();
2473          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2474          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2475            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2476              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2477              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2478              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2479              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2480              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2481              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2482    
2483            }
2484          }
2485    
2486        }
2487        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2488    
2489          // Borrow DataTagged input from Data object
2490          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2491    
2492          // Prepare the DataConstant input
2493          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2494    
2495          // Prepare a DataTagged output 2
2496          res = Data(0.0, shape1, arg_0_Z.getFunctionSpace());      // DataTagged output
2497          res.tag();
2498          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2499    
2500          // Prepare offset into DataConstant
2501          int offset_1 = tmp_1->getPointOffset(0,0);
2502          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2503    
2504          // Get the pointers to the actual data
2505          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2506          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2507    
2508    
2509          // Compute a result for the default
2510          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2511          // Compute a result for each tag
2512          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2513          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2514          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2515            tmp_2->addTag(i->first);
2516            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2517            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2518    
2519            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2520          }
2521    
2522        }
2523        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2524    
2525          // Borrow DataTagged input from Data object
2526          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2527    
2528          // Borrow DataTagged input from Data object
2529          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2530    
2531          // Prepare a DataTagged output 2
2532          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
2533          res.tag();        // DataTagged output
2534          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2535    
2536          // Get the pointers to the actual data
2537          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2538          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2539          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2540    
2541          // Compute a result for the default
2542          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2543          // Merge the tags
2544          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2545          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2546          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2547          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2548            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2549          }
2550          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2551            tmp_2->addTag(i->first);
2552          }
2553          // Compute a result for each tag
2554          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2555          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2556            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2557            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2558            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2559    
2560            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2561          }
2562    
2563        }
2564        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2565    
2566          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2567          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2568          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2569          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2570          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2571    
2572          int sampleNo_0,dataPointNo_0;
2573          int numSamples_0 = arg_0_Z.getNumSamples();
2574          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2575          res.requireWrite();
2576          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2577          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2578            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2579            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2580            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2581              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2582              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2583              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2584              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2585              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2586            }
2587          }
2588    
2589        }
2590        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2591          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2592          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2593          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2594          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2595    
2596          int sampleNo_0,dataPointNo_0;
2597          int numSamples_0 = arg_0_Z.getNumSamples();
2598          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2599          int offset_1 = tmp_1->getPointOffset(0,0);
2600          res.requireWrite();
2601          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2602          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2603            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2604              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2605              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2606              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2607              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2608              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2609              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2610            }
2611          }
2612    
2613    
2614        }
2615        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2616    
2617          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2618          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2619          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2620          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2621          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2622    
2623          int sampleNo_0,dataPointNo_0;
2624          int numSamples_0 = arg_0_Z.getNumSamples();
2625          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2626          res.requireWrite();
2627          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2628          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2629            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2630            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2631            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2632              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2633              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2634              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2635              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2636              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2637            }
2638          }
2639    
2640        }
2641        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2642    
2643          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2644          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2645          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2646          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2647          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2648    
2649          int sampleNo_0,dataPointNo_0;
2650          int numSamples_0 = arg_0_Z.getNumSamples();
2651          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2652          res.requireWrite();
2653          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2654          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2655            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2656              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2657              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2658              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2659              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2660              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2661              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2662              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2663            }
2664          }
2665    
2666        }
2667        else {
2668          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2669        }
2670    
2671      } else if (0 == rank1) {
2672        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2673          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2674          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2675          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2676          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2677          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2678        }
2679        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2680    
2681          // Prepare the DataConstant input
2682          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2683    
2684          // Borrow DataTagged input from Data object
2685          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2686    
2687          // Prepare a DataTagged output 2
2688          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2689          res.tag();
2690          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2691    
2692          // Prepare offset into DataConstant
2693          int offset_0 = tmp_0->getPointOffset(0,0);
2694          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2695    
2696          //Get the pointers to the actual data
2697          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2698          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2699    
2700          // Compute a result for the default
2701          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2702          // Compute a result for each tag
2703          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2704          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2705          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2706            tmp_2->addTag(i->first);
2707            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2708            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2709            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2710          }
2711        }
2712        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2713    
2714          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2715          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2716          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2717          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2718    
2719          int sampleNo_1,dataPointNo_1;
2720          int numSamples_1 = arg_1_Z.getNumSamples();
2721          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2722          int offset_0 = tmp_0->getPointOffset(0,0);
2723          res.requireWrite();
2724          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2725          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2726            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2727              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2728              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2729              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2730              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2731              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2732              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2733            }
2734          }
2735    
2736        }
2737        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2738    
2739          // Borrow DataTagged input from Data object
2740          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2741    
2742          // Prepare the DataConstant input
2743          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2744    
2745          // Prepare a DataTagged output 2
2746          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2747          res.tag();
2748          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2749    
2750          // Prepare offset into DataConstant
2751          int offset_1 = tmp_1->getPointOffset(0,0);
2752          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2753          // Get the pointers to the actual data
2754          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2755          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2756          // Compute a result for the default
2757          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2758          // Compute a result for each tag
2759          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2760          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2761          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2762            tmp_2->addTag(i->first);
2763            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2764            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2765            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2766          }
2767    
2768        }
2769        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2770    
2771          // Borrow DataTagged input from Data object
2772          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2773    
2774          // Borrow DataTagged input from Data object
2775          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2776    
2777          // Prepare a DataTagged output 2
2778          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2779          res.tag();        // DataTagged output
2780          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2781    
2782          // Get the pointers to the actual data
2783          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2784          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2785          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2786    
2787          // Compute a result for the default
2788          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2789          // Merge the tags
2790          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2791          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2792          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2793          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2794            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2795          }
2796          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2797            tmp_2->addTag(i->first);
2798          }
2799          // Compute a result for each tag
2800          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2801          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2802            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2803            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2804            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2805            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2806          }
2807    
2808        }
2809        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2810    
2811          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2812          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2813          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2814          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2815          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2816    
2817          int sampleNo_0,dataPointNo_0;
2818          int numSamples_0 = arg_0_Z.getNumSamples();
2819          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2820          res.requireWrite();
2821          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2822          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2823            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2824            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2825            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2826              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2827              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2828              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2829              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2830              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2831            }
2832          }
2833    
2834        }
2835        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2836          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2837          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2838          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2839          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2840    
2841          int sampleNo_0,dataPointNo_0;
2842          int numSamples_0 = arg_0_Z.getNumSamples();
2843          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2844          int offset_1 = tmp_1->getPointOffset(0,0);
2845          res.requireWrite();
2846          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2847          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2848            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2849              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2850              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2851              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2852              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2853              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2854              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2855            }
2856          }
2857    
2858    
2859        }
2860        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2861    
2862          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2863          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2864          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2865          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2866          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2867    
2868          int sampleNo_0,dataPointNo_0;
2869          int numSamples_0 = arg_0_Z.getNumSamples();
2870          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2871          res.requireWrite();
2872          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2873          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2874            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2875            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2876            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2877              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2878              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2879              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2880              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2881              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2882            }
2883          }
2884    
2885        }
2886        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2887    
2888          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2889          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2890          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2891          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2892          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2893    
2894          int sampleNo_0,dataPointNo_0;
2895          int numSamples_0 = arg_0_Z.getNumSamples();
2896          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2897          res.requireWrite();
2898          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2899          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2900            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2901              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2902              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2903              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2904              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2905              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2906              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2907              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2908            }
2909          }
2910    
2911        }
2912        else {
2913          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2914        }
2915    
2916      } else {
2917        throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible shapes");
2918      }
2919    
2920      return res;
2921    }
2922    
2923    template <typename UnaryFunction>
2924    Data
2925    C_TensorUnaryOperation(Data const &arg_0,
2926                           UnaryFunction operation)
2927    {
2928      if (arg_0.isEmpty())  // do this before we attempt to interpolate
2929      {
2930         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2931      }
2932      if (arg_0.isLazy())
2933      {
2934         throw DataException("Error - Operations not permitted on lazy data.");
2935      }
2936      // Interpolate if necessary and find an appropriate function space
2937      Data arg_0_Z = Data(arg_0);
2938    
2939      // Get rank and shape of inputs
2940      const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2941      int size0 = arg_0_Z.getDataPointSize();
2942    
2943      // Declare output Data object
2944      Data res;
2945    
2946      if (arg_0_Z.isConstant()) {
2947        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
2948        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2949        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2950        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2951      }
2952      else if (arg_0_Z.isTagged()) {
2953    
2954        // Borrow DataTagged input from Data object
2955        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2956    
2957        // Prepare a DataTagged output 2
2958        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());   // DataTagged output
2959        res.tag();
2960        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2961    
2962        // Get the pointers to the actual data
2963        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2964        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2965        // Compute a result for the default
2966        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2967        // Compute a result for each tag
2968        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2969        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2970        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2971          tmp_2->addTag(i->first);
2972          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2973          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2974          tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2975        }
2976    
2977      }
2978      else if (arg_0_Z.isExpanded()) {
2979    
2980        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
2981        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2982        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2983    
2984        int sampleNo_0,dataPointNo_0;
2985        int numSamples_0 = arg_0_Z.getNumSamples();
2986        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2987        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2988        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2989          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2990            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2991            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2992            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2993            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2994            tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2995          }
2996        }
2997      }
2998      else {
2999        throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
3000      }
3001    
3002      return res;
3003  }  }
3004    
3005  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26