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

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

  ViewVC Help
Powered by ViewVC 1.1.26