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

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

  ViewVC Help
Powered by ViewVC 1.1.26