/[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 2771 by jfenwick, Wed Nov 25 01:38:49 2009 UTC
# Line 1  Line 1 
 // $Id$  
 /*=============================================================================  
1    
2   ******************************************************************************  /*******************************************************
3   *                                                                            *  *
4   *       COPYRIGHT ACcESS 2004 -  All Rights Reserved                         *  * Copyright (c) 2003-2009 by University of Queensland
5   *                                                                            *  * Earth Systems Science Computational Center (ESSCC)
6   * This software is the property of ACcESS.  No part of this code             *  * http://www.uq.edu.au/esscc
7   * may be copied in any form or by any means without the expressed written    *  *
8   * consent of ACcESS.  Copying, use or modification of this software          *  * Primary Business: Queensland, Australia
9   * by any unauthorised person is illegal unless that                          *  * Licensed under the Open Software License version 3.0
10   * person has a software license agreement with ACcESS.                       *  * http://www.opensource.org/licenses/osl-3.0.php
11   *                                                                            *  *
12   ******************************************************************************  *******************************************************/
13    
14  ******************************************************************************/  
15    /** \file Data.h */
16    
17  #ifndef DATA_H  #ifndef DATA_H
18  #define DATA_H  #define DATA_H
19    #include "system_dep.h"
20    
21    #include "DataTypes.h"
22    #include "DataAbstract.h"
23    #include "DataAlgorithm.h"
24    #include "FunctionSpace.h"
25    #include "BinaryOp.h"
26    #include "UnaryOp.h"
27    #include "DataException.h"
28    
 #include "escript/Data/DataAbstract.h"  
 #include "escript/Data/DataTagged.h"  
 #include "escript/Data/FunctionSpace.h"  
 #include "escript/Data/BinaryOp.h"  
 #include "escript/Data/UnaryOp.h"  
29    
30  extern "C" {  extern "C" {
31  #include "escript/Data/DataC.h"  #include "DataC.h"
32    //#include <omp.h>
33  }  }
34    
35  #include <iostream>  #ifdef _OPENMP
36    #include <omp.h>
37    #endif
38    
39    #include "esysmpi.h"
40  #include <string>  #include <string>
 #include <memory>  
41  #include <algorithm>  #include <algorithm>
42    #include <sstream>
43    
44  #include <boost/shared_ptr.hpp>  #include <boost/shared_ptr.hpp>
45  #include <boost/python/object.hpp>  #include <boost/python/object.hpp>
 #include <boost/python/list.hpp>  
46  #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.  
 */  
47    
48  namespace escript {  namespace escript {
49    
50  //  //
51  // Forward declaration for various implimentations of Data.  // Forward declaration for various implementations of Data.
 class DataEmpty;  
52  class DataConstant;  class DataConstant;
53  class DataTagged;  class DataTagged;
54  class DataExpanded;  class DataExpanded;
55    class DataLazy;
56    
57    /**
58       \brief
59       Data represents a collection of datapoints.
60    
61       Description:
62       Internally, the datapoints are actually stored by a DataAbstract object.
63       The specific instance of DataAbstract used may vary over the lifetime
64       of the Data object.
65       Some methods on this class return references (eg getShape()).
66       These references should not be used after an operation which changes the underlying DataAbstract object.
67       Doing so will lead to invalid memory access.
68       This should not affect any methods exposed via boost::python.
69    */
70  class Data {  class Data {
71    
72    public:    public:
73    
74      // These typedefs allow function names to be cast to pointers
75      // to functions of the appropriate type when calling unaryOp etc.
76    typedef double (*UnaryDFunPtr)(double);    typedef double (*UnaryDFunPtr)(double);
77    typedef double (*BinaryDFunPtr)(double,double);    typedef double (*BinaryDFunPtr)(double,double);
78    
79    
80    /**    /**
81       Constructors.       Constructors.
82    */    */
# Line 80  class Data { Line 86  class Data {
86       Default constructor.       Default constructor.
87       Creates a DataEmpty object.       Creates a DataEmpty object.
88    */    */
89      ESCRIPT_DLL_API
90    Data();    Data();
91    
92    /**    /**
# Line 87  class Data { Line 94  class Data {
94       Copy constructor.       Copy constructor.
95       WARNING: Only performs a shallow copy.       WARNING: Only performs a shallow copy.
96    */    */
97      ESCRIPT_DLL_API
98    Data(const Data& inData);    Data(const Data& inData);
99    
100    /**    /**
# Line 95  class Data { Line 103  class Data {
103       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,
104       otherwise a shallow copy of inData is returned.       otherwise a shallow copy of inData is returned.
105    */    */
106      ESCRIPT_DLL_API
107    Data(const Data& inData,    Data(const Data& inData,
108         const FunctionSpace& what);         const FunctionSpace& what);
109    
110    /**    /**
111       \brief      \brief Copy Data from an existing vector
112       Constructor which copies data from a DataArrayView.    */
113    
114       \param value - Input - Data value for a single point.    ESCRIPT_DLL_API
115       \param what - Input - A description of what this data represents.    Data(const DataTypes::ValueType& value,
116       \param expanded - Input - Flag, if true fill the entire container with           const DataTypes::ShapeType& shape,
117                         the value. Otherwise a more efficient storage                   const FunctionSpace& what=FunctionSpace(),
118                         mechanism will be used.                   bool expanded=false);
   */  
   Data(const DataArrayView& value,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
119    
120    /**    /**
121       \brief       \brief
122       Constructor which creates a Data from a DataArrayView shape.       Constructor which creates a Data with points having the specified shape.
123    
124       \param value - Input - Single value applied to all Data.       \param value - Input - Single value applied to all Data.
125       \param dataPointShape - Input - The shape of each data point.       \param dataPointShape - Input - The shape of each data point.
# Line 123  class Data { Line 128  class Data {
128                         the given value. Otherwise a more efficient storage                         the given value. Otherwise a more efficient storage
129                         mechanism will be used.                         mechanism will be used.
130    */    */
131      ESCRIPT_DLL_API
132    Data(double value,    Data(double value,
133         const DataArrayView::ShapeType& dataPointShape=DataArrayView::ShapeType(),         const DataTypes::ShapeType& dataPointShape=DataTypes::ShapeType(),
134         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
135         bool expanded=false);         bool expanded=false);
136    
# Line 135  class Data { Line 141  class Data {
141       \param inData - Input - Input Data object.       \param inData - Input - Input Data object.
142       \param region - Input - Region to copy.       \param region - Input - Region to copy.
143    */    */
144      ESCRIPT_DLL_API
145    Data(const Data& inData,    Data(const Data& inData,
146         const DataArrayView::RegionType& region);         const DataTypes::RegionType& region);
   
   /**  
      \brief  
      Constructor which will create Tagged data if expanded is false.  
      No attempt is made to ensure the tag keys match the tag keys  
      within the function space.  
   
      \param tagKeys - Input - List of tag values.  
      \param values - Input - List of values, one for each tag.  
      \param defaultValue - Input - A default value, used if tag doesn't exist.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the appropriate values.  
   */  
   Data(const DataTagged::TagListType& tagKeys,  
        const DataTagged::ValueListType& values,  
        const DataArrayView& defaultValue,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
   
   /**  
      \brief  
      Constructor which copies data from a python numarray.  
   
      \param value - Input - Data value for a single point.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the value. Otherwise a more efficient storage  
                        mechanism will be used.  
   */  
   Data(const boost::python::numeric::array& value,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
147    
148    /**    /**
149       \brief       \brief
150       Constructor which copies data from any object that can be converted into       Constructor which copies data from any object that can be treated like a python array/sequence.
      a python numarray.  
151    
152       \param value - Input - Input data.       \param value - Input - Input data.
153       \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 155  class Data {
155                         the value. Otherwise a more efficient storage                         the value. Otherwise a more efficient storage
156                         mechanism will be used.                         mechanism will be used.
157    */    */
158      ESCRIPT_DLL_API
159    Data(const boost::python::object& value,    Data(const boost::python::object& value,
160         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
161         bool expanded=false);         bool expanded=false);
# Line 189  class Data { Line 163  class Data {
163    /**    /**
164       \brief       \brief
165       Constructor which creates a DataConstant.       Constructor which creates a DataConstant.
166       Copies data from any object that can be converted       Copies data from any object that can be treated like a python array/sequence.
167       into a numarray. All other parameters are copied from other.       All other parameters are copied from other.
168    
169       \param value - Input - Input data.       \param value - Input - Input data.
170       \param other - Input - contains all other parameters.       \param other - Input - contains all other parameters.
171    */    */
172      ESCRIPT_DLL_API
173    Data(const boost::python::object& value,    Data(const boost::python::object& value,
174         const Data& other);         const Data& other);
175    
# Line 202  class Data { Line 177  class Data {
177       \brief       \brief
178       Constructor which creates a DataConstant of "shape" with constant value.       Constructor which creates a DataConstant of "shape" with constant value.
179    */    */
180    Data(double value,    ESCRIPT_DLL_API
181         const boost::python::tuple& shape=boost::python::make_tuple(),    Data(double value,
182           const boost::python::tuple& shape=boost::python::make_tuple(),
183         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
184         bool expanded=false);         bool expanded=false);
185    
186    
187    
188      /**
189        \brief Create a Data using an existing DataAbstract. Warning: The new object assumes ownership of the pointer!
190        Once you have passed the pointer, do not delete it.
191      */
192      ESCRIPT_DLL_API
193      explicit Data(DataAbstract* underlyingdata);
194    
195      /**
196        \brief Create a Data based on the supplied DataAbstract
197      */
198      ESCRIPT_DLL_API
199      explicit Data(DataAbstract_ptr underlyingdata);
200    
201    /**    /**
202       \brief       \brief
203       Perform a deep copy.       Destructor
204    */    */
205      ESCRIPT_DLL_API
206      ~Data();
207    
208      /**
209         \brief Make this object a deep copy of "other".
210      */
211      ESCRIPT_DLL_API
212    void    void
213    copy(const Data& other);    copy(const Data& other);
214    
215    /**    /**
216         \brief Return a pointer to a deep copy of this object.
217      */
218      ESCRIPT_DLL_API
219      Data
220      copySelf();
221    
222    
223      /**
224         \brief produce a delayed evaluation version of this Data.
225      */
226      ESCRIPT_DLL_API
227      Data
228      delay();
229    
230      /**
231         \brief convert the current data into lazy data.
232      */
233      ESCRIPT_DLL_API
234      void
235      delaySelf();
236    
237    
238      /**
239       Member access methods.       Member access methods.
240    */    */
241    
242    /**    /**
243       \brief       \brief
244         switches on update protection
245    
246      */
247      ESCRIPT_DLL_API
248      void
249      setProtection();
250    
251      /**
252         \brief
253         Returns true, if the data object is protected against update
254    
255      */
256      ESCRIPT_DLL_API
257      bool
258      isProtected() const;
259    
260    
261    /**
262       \brief
263       Return the value of a data point as a python tuple.
264    */
265      ESCRIPT_DLL_API
266      const boost::python::object
267      getValueOfDataPointAsTuple(int dataPointNo);
268    
269      /**
270         \brief
271         sets the values of a data-point from a python object on this process
272      */
273      ESCRIPT_DLL_API
274      void
275      setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object);
276    
277      /**
278         \brief
279         sets the values of a data-point from a array-like object on this process
280      */
281      ESCRIPT_DLL_API
282      void
283      setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);
284    
285      /**
286         \brief
287         sets the values of a data-point on this process
288      */
289      ESCRIPT_DLL_API
290      void
291      setValueOfDataPoint(int dataPointNo, const double);
292    
293      /**
294         \brief Return a data point across all processors as a python tuple.
295      */
296      ESCRIPT_DLL_API
297      const boost::python::object
298      getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo);
299    
300      /**
301         \brief
302         Return the tag number associated with the given data-point.
303    
304      */
305      ESCRIPT_DLL_API
306      int
307      getTagNumber(int dpno);
308    
309      /**
310         \brief
311       Return the C wrapper for the Data object.       Return the C wrapper for the Data object.
312    */    */
313      ESCRIPT_DLL_API
314    escriptDataC    escriptDataC
315    getDataC();    getDataC();
316    
317    
318    
319    /**    /**
320       \brief       \brief
321       Return the C wrapper for the Data object - const version.       Return the C wrapper for the Data object - const version.
322    */    */
323      ESCRIPT_DLL_API
324    escriptDataC    escriptDataC
325    getDataC() const;    getDataC() const;
326    
   /**  
      \brief  
      Write the data as a string.  
   */  
   inline  
   std::string  
   toString() const  
   {  
     return m_data->toString();  
   }  
327    
328    /**    /**
329       \brief       \brief
330       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.  
331    */    */
332    inline    ESCRIPT_DLL_API
333    const DataArrayView&    std::string
334    getPointDataView() const    toString() const;
   {  
      return m_data->getPointDataView();  
   }  
335    
336    /**    /**
337       \brief       \brief
338       Whatever the current Data type make this into a DataExpanded.       Whatever the current Data type make this into a DataExpanded.
339    */    */
340      ESCRIPT_DLL_API
341    void    void
342    expand();    expand();
343    
# Line 269  class Data { Line 347  class Data {
347       Constant data to be converted to tagged. An attempt to convert       Constant data to be converted to tagged. An attempt to convert
348       Expanded data to tagged will throw an exception.       Expanded data to tagged will throw an exception.
349    */    */
350      ESCRIPT_DLL_API
351    void    void
352    tag();    tag();
353    
354    /**    /**
355        \brief If this data is lazy, then convert it to ready data.
356        What type of ready data depends on the expression. For example, Constant+Tagged==Tagged.
357      */
358      ESCRIPT_DLL_API
359      void
360      resolve();
361    
362    
363      /**
364       \brief Ensures data is ready for write access.
365      This means that the data will be resolved if lazy and will be copied if shared with another Data object.
366      \warning This method should only be called in single threaded sections of code. (It modifies m_data).
367      Do not create any Data objects from this one between calling requireWrite and getSampleDataRW.
368      Doing so might introduce additional sharing.
369      */
370      ESCRIPT_DLL_API
371      void
372      requireWrite();
373    
374      /**
375       \brief       \brief
376       Return true if this Data is expanded.       Return true if this Data is expanded.
377         \note To determine if a sample will contain separate values for each datapoint. Use actsExpanded instead.
378    */    */
379      ESCRIPT_DLL_API
380    bool    bool
381    isExpanded() const;    isExpanded() const;
382    
383    /**    /**
384       \brief       \brief
385         Return true if this Data is expanded or resolves to expanded.
386         That is, if it has a separate value for each datapoint in the sample.
387      */
388      ESCRIPT_DLL_API
389      bool
390      actsExpanded() const;
391      
392    
393      /**
394         \brief
395       Return true if this Data is tagged.       Return true if this Data is tagged.
396    */    */
397      ESCRIPT_DLL_API
398    bool    bool
399    isTagged() const;    isTagged() const;
400    
# Line 290  class Data { Line 402  class Data {
402       \brief       \brief
403       Return true if this Data is constant.       Return true if this Data is constant.
404    */    */
405      ESCRIPT_DLL_API
406    bool    bool
407    isConstant() const;    isConstant() const;
408    
409    /**    /**
410         \brief Return true if this Data is lazy.
411      */
412      ESCRIPT_DLL_API
413      bool
414      isLazy() const;
415    
416      /**
417         \brief Return true if this data is ready.
418      */
419      ESCRIPT_DLL_API
420      bool
421      isReady() const;
422    
423      /**
424       \brief       \brief
425       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
426    contains datapoints.
427    */    */
428      ESCRIPT_DLL_API
429    bool    bool
430    isEmpty() const;    isEmpty() const;
431    
# Line 304  class Data { Line 433  class Data {
433       \brief       \brief
434       Return the function space.       Return the function space.
435    */    */
436      ESCRIPT_DLL_API
437    inline    inline
438    const FunctionSpace&    const FunctionSpace&
439    getFunctionSpace() const    getFunctionSpace() const
# Line 315  class Data { Line 445  class Data {
445       \brief       \brief
446       Return a copy of the function space.       Return a copy of the function space.
447    */    */
448      ESCRIPT_DLL_API
449    const FunctionSpace    const FunctionSpace
450    getCopyOfFunctionSpace() const;    getCopyOfFunctionSpace() const;
451    
# Line 322  class Data { Line 453  class Data {
453       \brief       \brief
454       Return the domain.       Return the domain.
455    */    */
456      ESCRIPT_DLL_API
457    inline    inline
458    const AbstractDomain&  //   const AbstractDomain&
459      const_Domain_ptr
460    getDomain() const    getDomain() const
461    {    {
462       return getFunctionSpace().getDomain();       return getFunctionSpace().getDomain();
463    }    }
464    
465    
466      /**
467         \brief
468         Return the domain.
469         TODO: For internal use only.   This should be removed.
470      */
471      ESCRIPT_DLL_API
472      inline
473    //   const AbstractDomain&
474      Domain_ptr
475      getDomainPython() const
476      {
477         return getFunctionSpace().getDomainPython();
478      }
479    
480    /**    /**
481       \brief       \brief
482       Return a copy of the domain.       Return a copy of the domain.
483    */    */
484      ESCRIPT_DLL_API
485    const AbstractDomain    const AbstractDomain
486    getCopyOfDomain() const;    getCopyOfDomain() const;
487    
# Line 340  class Data { Line 489  class Data {
489       \brief       \brief
490       Return the rank of the point data.       Return the rank of the point data.
491    */    */
492      ESCRIPT_DLL_API
493    inline    inline
494    int    unsigned int
495    getDataPointRank() const    getDataPointRank() const
496    {    {
497      return m_data->getPointDataView().getRank();      return m_data->getRank();
498    }    }
499    
500    /**    /**
501       \brief       \brief
502         Return the number of data points
503      */
504      ESCRIPT_DLL_API
505      inline
506      int
507      getNumDataPoints() const
508      {
509        return getNumSamples() * getNumDataPointsPerSample();
510      }
511      /**
512         \brief
513       Return the number of samples.       Return the number of samples.
514    */    */
515      ESCRIPT_DLL_API
516    inline    inline
517    int    int
518    getNumSamples() const    getNumSamples() const
# Line 362  class Data { Line 524  class Data {
524       \brief       \brief
525       Return the number of data points per sample.       Return the number of data points per sample.
526    */    */
527      ESCRIPT_DLL_API
528    inline    inline
529    int    int
530    getNumDataPointsPerSample() const    getNumDataPointsPerSample() const
# Line 369  class Data { Line 532  class Data {
532      return m_data->getNumDPPSample();      return m_data->getNumDPPSample();
533    }    }
534    
535    
536      /**
537        \brief
538        Return the number of values in the shape for this object.
539      */
540      ESCRIPT_DLL_API
541      int
542      getNoValues() const
543      {
544        return m_data->getNoValues();
545      }
546    
547    
548      /**
549         \brief
550         dumps the object into a netCDF file
551      */
552      ESCRIPT_DLL_API
553      void
554      dump(const std::string fileName) const;
555    
556     /**
557      \brief returns the values of the object as a list of tuples (one for each datapoint).
558    
559      \param scalarastuple If true, scalar data will produce single valued tuples [(1,) (2,) ...]
560    If false, the result is a list of scalars [1, 2, ...]
561     */
562      ESCRIPT_DLL_API
563      const boost::python::object
564      toListOfTuples(bool scalarastuple=true);
565    
566    
567     /**
568        \brief
569        Return the sample data for the given sample no. This is not the
570        preferred interface but is provided for use by C code.
571        The bufferg parameter is only required for LazyData.
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.
611       NOTE: Construction of the DataArrayView is a relatively expensive       \param sampleNo - Input -
612       operation.       \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 -       \param sampleNo - Input -
622       \param dataPointNo - Input -       \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    
1760    };
1761    
1762    }   // end namespace escript
1763    
1764    
1765    // No, this is not supposed to be at the top of the file
1766    // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1767    // so that I can dynamic cast between them below.
1768    #include "DataReady.h"
1769    #include "DataLazy.h"
1770    
1771    namespace escript
1772    {
1773    
1774    inline
1775    const DataReady*
1776    Data::getReady() const
1777    {
1778       const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1779       EsysAssert((dr!=0), "Error - casting to DataReady.");
1780       return dr;
1781    }
1782    
1783    inline
1784    DataReady*
1785    Data::getReady()
1786    {
1787       DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1788       EsysAssert((dr!=0), "Error - casting to DataReady.");
1789       return dr;
1790    }
1791    
1792    // Be wary of using this for local operations since it (temporarily) increases reference count.
1793    // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1794    // getReady() instead
1795    inline
1796    DataReady_ptr
1797    Data::getReadyPtr()
1798    {
1799       DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1800       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1801       return dr;
1802    }
1803    
1804    
1805    inline
1806    const_DataReady_ptr
1807    Data::getReadyPtr() const
1808    {
1809       const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1810       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1811       return dr;
1812    }
1813    
1814    inline
1815    DataAbstract::ValueType::value_type*
1816    Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1817    {
1818       if (isLazy())
1819       {
1820        throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1821       }
1822       return getReady()->getSampleDataRW(sampleNo);
1823    }
1824    
1825    inline
1826    const DataAbstract::ValueType::value_type*
1827    Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo)
1828    {
1829       DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1830       if (l!=0)
1831       {
1832        size_t offset=0;
1833        const DataTypes::ValueType* res=l->resolveSample(sampleNo,offset);
1834        return &((*res)[offset]);
1835       }
1836       return getReady()->getSampleDataRO(sampleNo);
1837  }  }
1838    
1839    
1840    
1841    /**
1842       Modify a filename for MPI parallel output to multiple files
1843    */
1844    char *Escript_MPI_appendRankToFileName(const char *, int, int);
1845    
1846  /**  /**
1847     Binary Data object operators.     Binary Data object operators.
1848  */  */
1849    inline double rpow(double x,double y)
1850    {
1851        return pow(y,x);
1852    }
1853    
1854  /**  /**
1855    \brief    \brief
1856    Operator+    Operator+
1857    Takes two Data objects.    Takes two Data objects.
1858  */  */
1859  Data operator+(const Data& left, const Data& right);  ESCRIPT_DLL_API Data operator+(const Data& left, const Data& right);
1860    
1861  /**  /**
1862    \brief    \brief
1863    Operator-    Operator-
1864    Takes two Data objects.    Takes two Data objects.
1865  */  */
1866  Data operator-(const Data& left, const Data& right);  ESCRIPT_DLL_API Data operator-(const Data& left, const Data& right);
1867    
1868  /**  /**
1869    \brief    \brief
1870    Operator*    Operator*
1871    Takes two Data objects.    Takes two Data objects.
1872  */  */
1873  Data operator*(const Data& left, const Data& right);  ESCRIPT_DLL_API Data operator*(const Data& left, const Data& right);
1874    
1875  /**  /**
1876    \brief    \brief
1877    Operator/    Operator/
1878    Takes two Data objects.    Takes two Data objects.
1879  */  */
1880  Data operator/(const Data& left, const Data& right);  ESCRIPT_DLL_API Data operator/(const Data& left, const Data& right);
1881    
1882  /**  /**
1883    \brief    \brief
# Line 957  Data operator/(const Data& left, const D Line 1885  Data operator/(const Data& left, const D
1885    Takes LHS Data object and RHS python::object.    Takes LHS Data object and RHS python::object.
1886    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1887  */  */
1888  Data operator+(const Data& left, const boost::python::object& right);  ESCRIPT_DLL_API Data operator+(const Data& left, const boost::python::object& right);
1889    
1890  /**  /**
1891    \brief    \brief
# Line 965  Data operator+(const Data& left, const b Line 1893  Data operator+(const Data& left, const b
1893    Takes LHS Data object and RHS python::object.    Takes LHS Data object and RHS python::object.
1894    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1895  */  */
1896  Data operator-(const Data& left, const boost::python::object& right);  ESCRIPT_DLL_API Data operator-(const Data& left, const boost::python::object& right);
1897    
1898  /**  /**
1899    \brief    \brief
# Line 973  Data operator-(const Data& left, const b Line 1901  Data operator-(const Data& left, const b
1901    Takes LHS Data object and RHS python::object.    Takes LHS Data object and RHS python::object.
1902    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1903  */  */
1904  Data operator*(const Data& left, const boost::python::object& right);  ESCRIPT_DLL_API Data operator*(const Data& left, const boost::python::object& right);
1905    
1906  /**  /**
1907    \brief    \brief
# Line 981  Data operator*(const Data& left, const b Line 1909  Data operator*(const Data& left, const b
1909    Takes LHS Data object and RHS python::object.    Takes LHS Data object and RHS python::object.
1910    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1911  */  */
1912  Data operator/(const Data& left, const boost::python::object& right);  ESCRIPT_DLL_API Data operator/(const Data& left, const boost::python::object& right);
1913    
1914  /**  /**
1915    \brief    \brief
# Line 989  Data operator/(const Data& left, const b Line 1917  Data operator/(const Data& left, const b
1917    Takes LHS python::object and RHS Data object.    Takes LHS python::object and RHS Data object.
1918    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1919  */  */
1920  Data operator+(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator+(const boost::python::object& left, const Data& right);
1921    
1922  /**  /**
1923    \brief    \brief
# Line 997  Data operator+(const boost::python::obje Line 1925  Data operator+(const boost::python::obje
1925    Takes LHS python::object and RHS Data object.    Takes LHS python::object and RHS Data object.
1926    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1927  */  */
1928  Data operator-(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator-(const boost::python::object& left, const Data& right);
1929    
1930  /**  /**
1931    \brief    \brief
# Line 1005  Data operator-(const boost::python::obje Line 1933  Data operator-(const boost::python::obje
1933    Takes LHS python::object and RHS Data object.    Takes LHS python::object and RHS Data object.
1934    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1935  */  */
1936  Data operator*(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator*(const boost::python::object& left, const Data& right);
1937    
1938  /**  /**
1939    \brief    \brief
# Line 1013  Data operator*(const boost::python::obje Line 1941  Data operator*(const boost::python::obje
1941    Takes LHS python::object and RHS Data object.    Takes LHS python::object and RHS Data object.
1942    python::object must be convertable to Data type.    python::object must be convertable to Data type.
1943  */  */
1944  Data operator/(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator/(const boost::python::object& left, const Data& right);
1945    
1946    
1947    
1948  /**  /**
1949    \brief    \brief
1950    Output operator    Output operator
1951  */  */
1952  std::ostream& operator<<(std::ostream& o, const Data& data);  ESCRIPT_DLL_API std::ostream& operator<<(std::ostream& o, const Data& data);
1953    
1954  /**  /**
1955    \brief    \brief
1956    Return true if operands are equivalent, else return false.    Compute a tensor product of two Data objects
1957    NB: this operator does very little at this point, and isn't to    \param arg_0 - Input - Data object
1958    be relied on. Requires further implementation.    \param arg_1 - Input - Data object
1959      \param axis_offset - Input - axis offset
1960      \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1961  */  */
1962  //bool operator==(const Data& left, const Data& right);  ESCRIPT_DLL_API
1963    Data
1964    C_GeneralTensorProduct(Data& arg_0,
1965                         Data& arg_1,
1966                         int axis_offset=0,
1967                         int transpose=0);
1968    
1969  /**  /**
1970    \brief    \brief
# Line 1042  Data::binaryOp(const Data& right, Line 1979  Data::binaryOp(const Data& right,
1979  {  {
1980     //     //
1981     // 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
1982     if (getPointDataView().getRank()==0 && right.getPointDataView().getRank()!=0) {     if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1983       reshapeDataPoint(right.getPointDataView().getShape());       throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
1984       }
1985    
1986       if (isLazy() || right.isLazy())
1987       {
1988         throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
1989     }     }
1990     //     //
1991     // initially make the temporary a shallow copy     // initially make the temporary a shallow copy
1992     Data tempRight(right);     Data tempRight(right);
1993    
1994     if (getFunctionSpace()!=right.getFunctionSpace()) {     if (getFunctionSpace()!=right.getFunctionSpace()) {
1995       if (right.probeInterpolation(getFunctionSpace())) {       if (right.probeInterpolation(getFunctionSpace())) {
1996         //         //
1997         // an interpolation is required so create a new Data         // an interpolation is required so create a new Data
1998         tempRight=Data(right,this->getFunctionSpace());         tempRight=Data(right,this->getFunctionSpace());
1999       } else if (probeInterpolation(right.getFunctionSpace())) {       } else if (probeInterpolation(right.getFunctionSpace())) {
2000         //         //
2001         // interpolate onto the RHS function space         // interpolate onto the RHS function space
2002         Data tempLeft(*this,right.getFunctionSpace());         Data tempLeft(*this,right.getFunctionSpace());
2003         m_data=tempLeft.m_data;  //        m_data=tempLeft.m_data;
2004           set_m_data(tempLeft.m_data);
2005       }       }
2006     }     }
2007     operandCheck(tempRight);     operandCheck(tempRight);
# Line 1073  Data::binaryOp(const Data& right, Line 2017  Data::binaryOp(const Data& right,
2017       // of any data type       // of any data type
2018       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2019       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2020       escript::binaryOp(*leftC,*(tempRight.m_data.get()),operation);       escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
2021     } else if (isTagged()) {     } else if (isTagged()) {
2022       //       //
2023       // 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 2043  Data::binaryOp(const Data& right,
2043    
2044  /**  /**
2045    \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  
2046    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.
2047    Given operation combines each element of each data point, thus argument    Given operation combines each element of each data point, thus argument
2048    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.
2049    Calls escript::algorithm.    Calls escript::algorithm.
2050  */  */
2051  template <class UnaryFunction>  template <class BinaryFunction>
2052  inline  inline
2053  double  double
2054  Data::algorithm(UnaryFunction operation) const  Data::algorithm(BinaryFunction operation, double initial_value) const
2055  {  {
2056    if (isExpanded()) {    if (isExpanded()) {
2057      DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2058      EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");      EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2059      return escript::algorithm(*leftC,operation);      return escript::algorithm(*leftC,operation,initial_value);
2060    } else if (isTagged()) {    } else if (isTagged()) {
2061      DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
2062      EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");      EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");
2063      return escript::algorithm(*leftC,operation);      return escript::algorithm(*leftC,operation,initial_value);
2064    } else if (isConstant()) {    } else if (isConstant()) {
2065      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2066      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2067      return escript::algorithm(*leftC,operation);      return escript::algorithm(*leftC,operation,initial_value);
2068      } else if (isEmpty()) {
2069        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2070      } else if (isLazy()) {
2071        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2072      } else {
2073        throw DataException("Error - Data encapsulates an unknown type.");
2074    }    }
   return 0;  
2075  }  }
2076    
2077  /**  /**
2078    \brief    \brief
2079    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.
2080    Given operation combines each element within each data point into a scalar,    Given operation combines each element within each data point into a scalar,
2081    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
2082    rank 0 Data object.    rank 0 Data object.
2083    Calls escript::dp_algorithm.    Calls escript::dp_algorithm.
2084  */  */
2085  template <class UnaryFunction>  template <class BinaryFunction>
2086  inline  inline
2087  Data  Data
2088  dp_algorithm(const Data& data,  Data::dp_algorithm(BinaryFunction operation, double initial_value) const
              UnaryFunction operation)  
2089  {  {
2090    Data result(0,DataArrayView::ShapeType(),data.getFunctionSpace(),data.isExpanded());    if (isEmpty()) {
2091    if (data.isExpanded()) {      throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2092      DataExpanded* dataE=dynamic_cast<DataExpanded*>(data.m_data.get());    }
2093      else if (isExpanded()) {
2094        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2095        DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2096      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2097      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2098      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2099      escript::dp_algorithm(*dataE,*resultE,operation);      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2100    } else if (data.isTagged()) {      return result;
2101      DataTagged* dataT=dynamic_cast<DataTagged*>(data.m_data.get());    }
2102      DataTagged* resultT=dynamic_cast<DataTagged*>(result.m_data.get());    else if (isTagged()) {
2103        DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());
2104      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2105      EsysAssert((resultT!=0), "Programming error - casting result to DataTagged.");      DataTypes::ValueType defval(1);
2106      escript::dp_algorithm(*dataT,*resultT,operation);      defval[0]=0;
2107    } else if (data.isConstant()) {      DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2108      DataConstant* dataC=dynamic_cast<DataConstant*>(data.m_data.get());      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2109        return Data(resultT);   // note: the Data object now owns the resultT pointer
2110      }
2111      else if (isConstant()) {
2112        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2113        DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2114      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2115      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2116      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2117      escript::dp_algorithm(*dataC,*resultC,operation);      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2118        return result;
2119      } else if (isLazy()) {
2120        throw DataException("Error - Operations not permitted on instances of DataLazy.");
2121      } else {
2122        throw DataException("Error - Data encapsulates an unknown type.");
2123      }
2124    }
2125    
2126    /**
2127      \brief
2128      Compute a tensor operation with two Data objects
2129      \param arg_0 - Input - Data object
2130      \param arg_1 - Input - Data object
2131      \param operation - Input - Binary op functor
2132    */
2133    template <typename BinaryFunction>
2134    inline
2135    Data
2136    C_TensorBinaryOperation(Data const &arg_0,
2137                            Data const &arg_1,
2138                            BinaryFunction operation)
2139    {
2140      if (arg_0.isEmpty() || arg_1.isEmpty())
2141      {
2142         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2143      }
2144      if (arg_0.isLazy() || arg_1.isLazy())
2145      {
2146         throw DataException("Error - Operations not permitted on lazy data.");
2147      }
2148      // Interpolate if necessary and find an appropriate function space
2149      Data arg_0_Z, arg_1_Z;
2150      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2151        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2152          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2153          arg_1_Z = Data(arg_1);
2154        }
2155        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2156          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2157          arg_0_Z =Data(arg_0);
2158        }
2159        else {
2160          throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible function spaces.");
2161        }
2162      } else {
2163          arg_0_Z = Data(arg_0);
2164          arg_1_Z = Data(arg_1);
2165      }
2166      // Get rank and shape of inputs
2167      int rank0 = arg_0_Z.getDataPointRank();
2168      int rank1 = arg_1_Z.getDataPointRank();
2169      DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2170      DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2171      int size0 = arg_0_Z.getDataPointSize();
2172      int size1 = arg_1_Z.getDataPointSize();
2173      // Declare output Data object
2174      Data res;
2175    
2176      if (shape0 == shape1) {
2177        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2178          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2179          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2180          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2181          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2182    
2183          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2184        }
2185        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2186    
2187          // Prepare the DataConstant input
2188          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2189    
2190          // Borrow DataTagged input from Data object
2191          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2192    
2193          // Prepare a DataTagged output 2
2194          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2195          res.tag();
2196          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2197    
2198          // Prepare offset into DataConstant
2199          int offset_0 = tmp_0->getPointOffset(0,0);
2200          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2201    
2202          // Get the pointers to the actual data
2203          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2204          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2205    
2206          // Compute a result for the default
2207          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2208          // Compute a result for each tag
2209          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2210          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2211          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2212            tmp_2->addTag(i->first);
2213            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2214            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2215    
2216            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2217          }
2218    
2219        }
2220        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2221          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2222          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2223          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2224          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2225    
2226          int sampleNo_1,dataPointNo_1;
2227          int numSamples_1 = arg_1_Z.getNumSamples();
2228          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2229          int offset_0 = tmp_0->getPointOffset(0,0);
2230          res.requireWrite();
2231          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2232          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2233            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2234              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2235              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2236              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2237              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2238              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2239              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2240            }
2241          }
2242    
2243        }
2244        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2245          // Borrow DataTagged input from Data object
2246          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2247    
2248          // Prepare the DataConstant input
2249          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2250    
2251          // Prepare a DataTagged output 2
2252          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2253          res.tag();
2254          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2255    
2256          // Prepare offset into DataConstant
2257          int offset_1 = tmp_1->getPointOffset(0,0);
2258    
2259          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2260          // Get the pointers to the actual data
2261          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2262          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2263          // Compute a result for the default
2264          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2265          // Compute a result for each tag
2266          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2267          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2268          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2269            tmp_2->addTag(i->first);
2270            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2271            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2272            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2273          }
2274    
2275        }
2276        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2277          // Borrow DataTagged input from Data object
2278          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2279    
2280          // Borrow DataTagged input from Data object
2281          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2282    
2283          // Prepare a DataTagged output 2
2284          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2285          res.tag();        // DataTagged output
2286          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2287    
2288          // Get the pointers to the actual data
2289          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2290          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2291          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2292    
2293          // Compute a result for the default
2294          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2295          // Merge the tags
2296          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2297          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2298          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2299          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2300            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2301          }
2302          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2303            tmp_2->addTag(i->first);
2304          }
2305          // Compute a result for each tag
2306          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2307          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2308    
2309            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2310            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2311            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2312    
2313            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2314          }
2315    
2316        }
2317        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2318          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2319          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2320          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2321          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2322          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2323    
2324          int sampleNo_0,dataPointNo_0;
2325          int numSamples_0 = arg_0_Z.getNumSamples();
2326          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2327          res.requireWrite();
2328          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2329          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2330            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2331            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2332            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2333              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2334              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2335              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2336              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2337              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2338            }
2339          }
2340    
2341        }
2342        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2343          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2344          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2345          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2346          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2347    
2348          int sampleNo_0,dataPointNo_0;
2349          int numSamples_0 = arg_0_Z.getNumSamples();
2350          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2351          int offset_1 = tmp_1->getPointOffset(0,0);
2352          res.requireWrite();
2353          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2354          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2355            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2356              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2357              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2358    
2359              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2360              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2361              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2362    
2363    
2364              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2365            }
2366          }
2367    
2368        }
2369        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2370          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2371          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2372          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2373          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2374          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2375    
2376          int sampleNo_0,dataPointNo_0;
2377          int numSamples_0 = arg_0_Z.getNumSamples();
2378          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2379          res.requireWrite();
2380          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2381          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2382            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2383            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2384            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2385              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2386              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2387              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2388              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2389              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2390            }
2391          }
2392    
2393        }
2394        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2395          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2396          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2397          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2398          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2399          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2400    
2401          int sampleNo_0,dataPointNo_0;
2402          int numSamples_0 = arg_0_Z.getNumSamples();
2403          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2404          res.requireWrite();
2405          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2406          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2407          dataPointNo_0=0;
2408    //        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2409              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2410              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2411              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2412              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2413              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2414              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2415              tensor_binary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_1, ptr_2, operation);
2416    //       }
2417          }
2418    
2419        }
2420        else {
2421          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2422        }
2423    
2424      } else if (0 == rank0) {
2425        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2426          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
2427          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2428          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2429          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2430          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2431        }
2432        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2433    
2434          // Prepare the DataConstant input
2435          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2436    
2437          // Borrow DataTagged input from Data object
2438          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2439    
2440          // Prepare a DataTagged output 2
2441          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataTagged output
2442          res.tag();
2443          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2444    
2445          // Prepare offset into DataConstant
2446          int offset_0 = tmp_0->getPointOffset(0,0);
2447          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2448    
2449          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2450          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2451    
2452          // Compute a result for the default
2453          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2454          // Compute a result for each tag
2455          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2456          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2457          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2458            tmp_2->addTag(i->first);
2459            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2460            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2461            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2462          }
2463    
2464        }
2465        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2466    
2467          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2468          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2469          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2470          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2471    
2472          int sampleNo_1,dataPointNo_1;
2473          int numSamples_1 = arg_1_Z.getNumSamples();
2474          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2475          int offset_0 = tmp_0->getPointOffset(0,0);
2476          res.requireWrite();
2477          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2478          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2479            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2480              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2481              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2482              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2483              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2484              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2485              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2486    
2487            }
2488          }
2489    
2490        }
2491        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2492    
2493          // Borrow DataTagged input from Data object
2494          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2495    
2496          // Prepare the DataConstant input
2497          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2498    
2499          // Prepare a DataTagged output 2
2500          res = Data(0.0, shape1, arg_0_Z.getFunctionSpace());      // DataTagged output
2501          res.tag();
2502          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2503    
2504          // Prepare offset into DataConstant
2505          int offset_1 = tmp_1->getPointOffset(0,0);
2506          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2507    
2508          // Get the pointers to the actual data
2509          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2510          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2511    
2512    
2513          // Compute a result for the default
2514          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2515          // Compute a result for each tag
2516          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2517          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2518          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2519            tmp_2->addTag(i->first);
2520            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2521            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2522    
2523            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2524          }
2525    
2526        }
2527        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2528    
2529          // Borrow DataTagged input from Data object
2530          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2531    
2532          // Borrow DataTagged input from Data object
2533          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2534    
2535          // Prepare a DataTagged output 2
2536          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
2537          res.tag();        // DataTagged output
2538          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2539    
2540          // Get the pointers to the actual data
2541          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2542          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2543          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2544    
2545          // Compute a result for the default
2546          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2547          // Merge the tags
2548          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2549          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2550          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2551          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2552            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2553          }
2554          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2555            tmp_2->addTag(i->first);
2556          }
2557          // Compute a result for each tag
2558          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2559          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2560            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2561            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2562            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2563    
2564            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2565          }
2566    
2567        }
2568        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2569    
2570          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2571          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2572          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2573          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2574          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2575    
2576          int sampleNo_0,dataPointNo_0;
2577          int numSamples_0 = arg_0_Z.getNumSamples();
2578          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2579          res.requireWrite();
2580          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2581          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2582            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2583            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2584            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2585              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2586              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2587              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2588              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2589              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2590            }
2591          }
2592    
2593        }
2594        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2595          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2596          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2597          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2598          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2599    
2600          int sampleNo_0,dataPointNo_0;
2601          int numSamples_0 = arg_0_Z.getNumSamples();
2602          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2603          int offset_1 = tmp_1->getPointOffset(0,0);
2604          res.requireWrite();
2605          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2606          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2607            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2608              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2609              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2610              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2611              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2612              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2613              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2614            }
2615          }
2616    
2617    
2618        }
2619        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2620    
2621          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2622          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2623          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2624          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2625          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2626    
2627          int sampleNo_0,dataPointNo_0;
2628          int numSamples_0 = arg_0_Z.getNumSamples();
2629          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2630          res.requireWrite();
2631          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2632          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2633            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2634            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2635            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2636              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2637              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2638              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2639              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2640              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2641            }
2642          }
2643    
2644        }
2645        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2646    
2647          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2648          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2649          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2650          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2651          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2652    
2653          int sampleNo_0,dataPointNo_0;
2654          int numSamples_0 = arg_0_Z.getNumSamples();
2655          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2656          res.requireWrite();
2657          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2658          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2659            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2660              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2661              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2662              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2663              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2664              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2665              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2666              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2667            }
2668          }
2669    
2670        }
2671        else {
2672          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2673        }
2674    
2675      } else if (0 == rank1) {
2676        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2677          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2678          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2679          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2680          double *ptr_2 = &(res.getDataAtOffsetRW(0));
2681          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2682        }
2683        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2684    
2685          // Prepare the DataConstant input
2686          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2687    
2688          // Borrow DataTagged input from Data object
2689          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2690    
2691          // Prepare a DataTagged output 2
2692          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2693          res.tag();
2694          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2695    
2696          // Prepare offset into DataConstant
2697          int offset_0 = tmp_0->getPointOffset(0,0);
2698          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2699    
2700          //Get the pointers to the actual data
2701          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2702          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2703    
2704          // Compute a result for the default
2705          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2706          // Compute a result for each tag
2707          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2708          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2709          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2710            tmp_2->addTag(i->first);
2711            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2712            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2713            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2714          }
2715        }
2716        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2717    
2718          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2719          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2720          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2721          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2722    
2723          int sampleNo_1,dataPointNo_1;
2724          int numSamples_1 = arg_1_Z.getNumSamples();
2725          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2726          int offset_0 = tmp_0->getPointOffset(0,0);
2727          res.requireWrite();
2728          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2729          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2730            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2731              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2732              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2733              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2734              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2735              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2736              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2737            }
2738          }
2739    
2740        }
2741        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2742    
2743          // Borrow DataTagged input from Data object
2744          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2745    
2746          // Prepare the DataConstant input
2747          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2748    
2749          // Prepare a DataTagged output 2
2750          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2751          res.tag();
2752          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2753    
2754          // Prepare offset into DataConstant
2755          int offset_1 = tmp_1->getPointOffset(0,0);
2756          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2757          // Get the pointers to the actual data
2758          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2759          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2760          // Compute a result for the default
2761          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2762          // Compute a result for each tag
2763          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2764          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2765          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2766            tmp_2->addTag(i->first);
2767            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2768            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2769            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2770          }
2771    
2772        }
2773        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2774    
2775          // Borrow DataTagged input from Data object
2776          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2777    
2778          // Borrow DataTagged input from Data object
2779          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2780    
2781          // Prepare a DataTagged output 2
2782          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2783          res.tag();        // DataTagged output
2784          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2785    
2786          // Get the pointers to the actual data
2787          const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2788          const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2789          double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2790    
2791          // Compute a result for the default
2792          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2793          // Merge the tags
2794          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2795          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2796          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2797          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2798            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2799          }
2800          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2801            tmp_2->addTag(i->first);
2802          }
2803          // Compute a result for each tag
2804          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2805          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2806            const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2807            const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2808            double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2809            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2810          }
2811    
2812        }
2813        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2814    
2815          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2816          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2817          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2818          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2819          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2820    
2821          int sampleNo_0,dataPointNo_0;
2822          int numSamples_0 = arg_0_Z.getNumSamples();
2823          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2824          res.requireWrite();
2825          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2826          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2827            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2828            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2829            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2830              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2831              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2832              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2833              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2834              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2835            }
2836          }
2837    
2838        }
2839        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2840          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2841          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2842          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2843          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2844    
2845          int sampleNo_0,dataPointNo_0;
2846          int numSamples_0 = arg_0_Z.getNumSamples();
2847          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2848          int offset_1 = tmp_1->getPointOffset(0,0);
2849          res.requireWrite();
2850          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2851          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2852            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2853              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2854              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2855              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2856              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2857              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2858              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2859            }
2860          }
2861    
2862    
2863        }
2864        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2865    
2866          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2867          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2868          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2869          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2870          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2871    
2872          int sampleNo_0,dataPointNo_0;
2873          int numSamples_0 = arg_0_Z.getNumSamples();
2874          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2875          res.requireWrite();
2876          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2877          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2878            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2879            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2880            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2881              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2882              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2883              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2884              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2885              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2886            }
2887          }
2888    
2889        }
2890        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2891    
2892          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2893          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2894          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2895          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2896          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2897    
2898          int sampleNo_0,dataPointNo_0;
2899          int numSamples_0 = arg_0_Z.getNumSamples();
2900          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2901          res.requireWrite();
2902          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2903          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2904            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2905              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2906              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2907              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2908              const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2909              const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2910              double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2911              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2912            }
2913          }
2914    
2915        }
2916        else {
2917          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2918        }
2919    
2920      } else {
2921        throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible shapes");
2922      }
2923    
2924      return res;
2925    }
2926    
2927    template <typename UnaryFunction>
2928    Data
2929    C_TensorUnaryOperation(Data const &arg_0,
2930                           UnaryFunction operation)
2931    {
2932      if (arg_0.isEmpty())  // do this before we attempt to interpolate
2933      {
2934         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2935      }
2936      if (arg_0.isLazy())
2937      {
2938         throw DataException("Error - Operations not permitted on lazy data.");
2939    }    }
2940    return result;    // Interpolate if necessary and find an appropriate function space
2941      Data arg_0_Z = Data(arg_0);
2942    
2943      // Get rank and shape of inputs
2944      const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2945      int size0 = arg_0_Z.getDataPointSize();
2946    
2947      // Declare output Data object
2948      Data res;
2949    
2950      if (arg_0_Z.isConstant()) {
2951        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
2952        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2953        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2954        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2955      }
2956      else if (arg_0_Z.isTagged()) {
2957    
2958        // Borrow DataTagged input from Data object
2959        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2960    
2961        // Prepare a DataTagged output 2
2962        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());   // DataTagged output
2963        res.tag();
2964        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2965    
2966        // Get the pointers to the actual data
2967        const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2968        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2969        // Compute a result for the default
2970        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2971        // Compute a result for each tag
2972        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2973        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2974        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2975          tmp_2->addTag(i->first);
2976          const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2977          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2978          tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2979        }
2980    
2981      }
2982      else if (arg_0_Z.isExpanded()) {
2983    
2984        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
2985        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2986        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2987    
2988        int sampleNo_0,dataPointNo_0;
2989        int numSamples_0 = arg_0_Z.getNumSamples();
2990        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2991        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2992        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2993        dataPointNo_0=0;
2994    //      for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2995            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2996            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2997            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2998            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2999            tensor_unary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_2, operation);
3000    //      }
3001        }
3002      }
3003      else {
3004        throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
3005      }
3006    
3007      return res;
3008  }  }
3009    
3010  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26