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

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

  ViewVC Help
Powered by ViewVC 1.1.26