/[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

revision 757 by woo409, Mon Jun 26 13:12:56 2006 UTC revision 2105 by jfenwick, Fri Nov 28 01:52:12 2008 UTC
# Line 1  Line 1 
1  // $Id$  
2  /*  /*******************************************************
3   ************************************************************  *
4   *          Copyright 2006 by ACcESS MNRF                   *  * Copyright (c) 2003-2008 by University of Queensland
5   *                                                          *  * Earth Systems Science Computational Center (ESSCC)
6   *              http://www.access.edu.au                    *  * http://www.uq.edu.au/esscc
7   *       Primary Business: Queensland, Australia            *  *
8   *  Licensed under the Open Software License version 3.0    *  * Primary Business: Queensland, Australia
9   *     http://www.opensource.org/licenses/osl-3.0.php       *  * Licensed under the Open Software License version 3.0
10   *                                                          *  * http://www.opensource.org/licenses/osl-3.0.php
11   ************************************************************  *
12  */  *******************************************************/
13    
14    
15  /** \file Data.h */  /** \file Data.h */
16    
# Line 23  Line 24 
24  #include "BinaryOp.h"  #include "BinaryOp.h"
25  #include "UnaryOp.h"  #include "UnaryOp.h"
26  #include "DataException.h"  #include "DataException.h"
27    #include "DataTypes.h"
28    
29  extern "C" {  extern "C" {
30  #include "DataC.h"  #include "DataC.h"
31  #include "paso/Paso.h"  /* #include "paso/Paso.h" doesn't belong in this file...causes trouble for BruceFactory.cpp */
32  }  }
33    
34    #include "esysmpi.h"
35  #include <string>  #include <string>
36  #include <algorithm>  #include <algorithm>
37    #include <sstream>
38    
39    
40  #include <boost/shared_ptr.hpp>  #include <boost/shared_ptr.hpp>
41  #include <boost/python/object.hpp>  #include <boost/python/object.hpp>
# Line 47  class DataExpanded; Line 52  class DataExpanded;
52    
53  /**  /**
54     \brief     \brief
55     Data creates the appropriate Data object for the given construction     Data represents a collection of datapoints.
    arguments.  
56    
57     Description:     Description:
58     Data is essentially a factory class which creates the appropriate Data     Internally, the datapoints are actually stored by a DataAbstract object.
59     object for the given construction arguments. It retains control over     The specific instance of DataAbstract used may vary over the lifetime
60     the object created for the lifetime of the object.     of the Data object.
61     The type of Data object referred to may change during the lifetime of     Some methods on this class return references (eg getShape()).
62     the Data object.     These references should not be used after an operation which changes the underlying DataAbstract object.
63       Doing so will lead to invalid memory access.
64       This should not affect any methods exposed via boost::python.
65  */  */
66  class Data {  class Data {
67    
# Line 66  class Data { Line 72  class Data {
72    typedef double (*UnaryDFunPtr)(double);    typedef double (*UnaryDFunPtr)(double);
73    typedef double (*BinaryDFunPtr)(double,double);    typedef double (*BinaryDFunPtr)(double,double);
74    
75    
76    /**    /**
77       Constructors.       Constructors.
78    */    */
# Line 97  class Data { Line 104  class Data {
104         const FunctionSpace& what);         const FunctionSpace& what);
105    
106    /**    /**
107       \brief      \brief Copy Data from an existing vector
108       Constructor which copies data from a DataArrayView.    */
109    
      \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.  
   */  
110    ESCRIPT_DLL_API    ESCRIPT_DLL_API
111    Data(const DataArrayView& value,    Data(const DataTypes::ValueType& value,
112         const FunctionSpace& what=FunctionSpace(),           const DataTypes::ShapeType& shape,
113         bool expanded=false);                   const FunctionSpace& what=FunctionSpace(),
114                     bool expanded=false);
115    
116    /**    /**
117       \brief       \brief
# Line 124  class Data { Line 126  class Data {
126    */    */
127    ESCRIPT_DLL_API    ESCRIPT_DLL_API
128    Data(double value,    Data(double value,
129         const DataArrayView::ShapeType& dataPointShape=DataArrayView::ShapeType(),         const DataTypes::ShapeType& dataPointShape=DataTypes::ShapeType(),
130         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
131         bool expanded=false);         bool expanded=false);
132    
# Line 137  class Data { Line 139  class Data {
139    */    */
140    ESCRIPT_DLL_API    ESCRIPT_DLL_API
141    Data(const Data& inData,    Data(const Data& inData,
142         const DataArrayView::RegionType& region);         const DataTypes::RegionType& region);
   
   /**  
      \brief  
      Constructor which will create Tagged data if expanded is false.  
      No attempt is made to ensure the tag keys match the tag keys  
      within the function space.  
   
      \param tagKeys - Input - List of tag values.  
      \param values - Input - List of values, one for each tag.  
      \param defaultValue - Input - A default value, used if tag doesn't exist.  
      \param what - Input - A description of what this data represents.  
      \param expanded - Input - Flag, if true fill the entire container with  
                        the appropriate values.  
     ==>*  
   */  
   ESCRIPT_DLL_API  
   Data(const DataTagged::TagListType& tagKeys,  
        const DataTagged::ValueListType& values,  
        const DataArrayView& defaultValue,  
        const FunctionSpace& what=FunctionSpace(),  
        bool expanded=false);  
143    
144    /**    /**
145       \brief       \brief
# Line 209  class Data { Line 190  class Data {
190       Constructor which creates a DataConstant of "shape" with constant value.       Constructor which creates a DataConstant of "shape" with constant value.
191    */    */
192    ESCRIPT_DLL_API    ESCRIPT_DLL_API
193    Data(double value,    Data(double value,
194         const boost::python::tuple& shape=boost::python::make_tuple(),         const boost::python::tuple& shape=boost::python::make_tuple(),
195         const FunctionSpace& what=FunctionSpace(),         const FunctionSpace& what=FunctionSpace(),
196         bool expanded=false);         bool expanded=false);
197    
198    
199    
200      /**
201        \brief Create a Data using an existing DataAbstract. Warning: The new object assumes ownership of the pointer!
202        Once you have passed the pointer, do not delete it.
203      */
204      ESCRIPT_DLL_API
205      explicit Data(DataAbstract* underlyingdata);
206    
207      /**
208        \brief Create a Data based on the supplied DataAbstract
209      */
210      ESCRIPT_DLL_API
211      explicit Data(DataAbstract_ptr underlyingdata);
212    
213    /**    /**
214       \brief       \brief
215       Destructor       Destructor
# Line 221  class Data { Line 218  class Data {
218    ~Data();    ~Data();
219    
220    /**    /**
221       \brief       \brief Make this object a deep copy of "other".
      Perform a deep copy.  
222    */    */
223    ESCRIPT_DLL_API    ESCRIPT_DLL_API
224    void    void
225    copy(const Data& other);    copy(const Data& other);
226    
227    /**    /**
228         \brief Return a pointer to a deep copy of this object.
229      */
230      ESCRIPT_DLL_API
231      Data
232      copySelf();
233    
234    
235      /**
236         \brief produce a delayed evaluation version of this Data.
237      */
238      ESCRIPT_DLL_API
239      Data
240      delay();
241    
242      /**
243         \brief convert the current data into lazy data.
244      */
245      ESCRIPT_DLL_API
246      void
247      delaySelf();
248    
249    
250      /**
251       Member access methods.       Member access methods.
252    */    */
253    
254    /**    /**
255       \brief       \brief
256       Return the values of all data-points as a single python numarray object.       switches on update protection
257    
258    */    */
259    ESCRIPT_DLL_API    ESCRIPT_DLL_API
260    const boost::python::numeric::array    void
261    convertToNumArray();    setProtection();
262    
263    /**    /**
264       \brief       \brief
265       Return the values of all data-points for the given sample as a single python numarray object.       Returns trueif the data object is protected against update
266    
267    */    */
268    ESCRIPT_DLL_API    ESCRIPT_DLL_API
269    const boost::python::numeric::array    bool
270    convertToNumArrayFromSampleNo(int sampleNo);    isProtected() const;
271    
272    /**    /**
273       \brief       \brief
274       Return the value of the specified data-point as a single python numarray object.       Return the values of a data point on this process
275    */    */
276    ESCRIPT_DLL_API    ESCRIPT_DLL_API
277    const boost::python::numeric::array    const boost::python::numeric::array
278    convertToNumArrayFromDPNo(int sampleNo,    getValueOfDataPoint(int dataPointNo);
279                              int dataPointNo);  
280      /**
281         \brief
282         sets the values of a data-point from a python object on this process
283      */
284      ESCRIPT_DLL_API
285      void
286      setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object);
287    
288      /**
289         \brief
290         sets the values of a data-point from a numarray object on this process
291      */
292      ESCRIPT_DLL_API
293      void
294      setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array&);
295    
296    /**    /**
297       \brief       \brief
298       Fills the expanded Data object from values of a python numarray object.       sets the values of a data-point on this process
299    */    */
300    ESCRIPT_DLL_API    ESCRIPT_DLL_API
301    void    void
302    fillFromNumArray(const boost::python::numeric::array);    setValueOfDataPoint(int dataPointNo, const double);
303    
304      /**
305         \brief
306         Return the value of the specified data-point across all processors
307      */
308      ESCRIPT_DLL_API
309      const boost::python::numeric::array
310      getValueOfGlobalDataPoint(int procNo, int dataPointNo);
311    
312    /**    /**
313       \brief       \brief
314       Return the tag number associated with the given data-point.       Return the tag number associated with the given data-point.
315    
      The data-point number here corresponds to the data-point number in the  
      numarray returned by convertToNumArray.  
316    */    */
317    ESCRIPT_DLL_API    ESCRIPT_DLL_API
318    int    int
# Line 284  class Data { Line 326  class Data {
326    escriptDataC    escriptDataC
327    getDataC();    getDataC();
328    
329    
330    
331    
332    
333    
334    // REMOVE ME
335    // ESCRIPT_DLL_API
336    // void
337    // CompareDebug(const Data& rd);
338    
339    
340    /**    /**
341       \brief       \brief
342       Return the C wrapper for the Data object - const version.       Return the C wrapper for the Data object - const version.
# Line 294  class Data { Line 347  class Data {
347    
348    /**    /**
349       \brief       \brief
350       Write the data as a string.       Write the data as a string. For large amounts of data, a summary is printed.
351    */    */
352    ESCRIPT_DLL_API    ESCRIPT_DLL_API
   inline  
353    std::string    std::string
354    toString() const    toString() const;
   {  
     return m_data->toString();  
   }  
   
   /**  
      \brief  
      Return the DataArrayView of the point data. This essentially contains  
      the shape information for each data point although it also may be used  
      to manipulate the point data.  
   */  
   ESCRIPT_DLL_API  
   inline  
   const DataArrayView&  
   getPointDataView() const  
   {  
      return m_data->getPointDataView();  
   }  
355    
356    /**    /**
357       \brief       \brief
# Line 331  class Data { Line 366  class Data {
366       If possible convert this Data to DataTagged. This will only allow       If possible convert this Data to DataTagged. This will only allow
367       Constant data to be converted to tagged. An attempt to convert       Constant data to be converted to tagged. An attempt to convert
368       Expanded data to tagged will throw an exception.       Expanded data to tagged will throw an exception.
     ==>*  
369    */    */
370    ESCRIPT_DLL_API    ESCRIPT_DLL_API
371    void    void
372    tag();    tag();
373    
374    /**    /**
375        \brief If this data is lazy, then convert it to ready data.
376        What type of ready data depends on the expression. For example, Constant+Tagged==Tagged.
377      */
378      ESCRIPT_DLL_API
379      void
380      resolve();
381    
382    
383      /**
384       \brief       \brief
385       Return true if this Data is expanded.       Return true if this Data is expanded.
386    */    */
# Line 362  class Data { Line 405  class Data {
405    isConstant() const;    isConstant() const;
406    
407    /**    /**
408         \brief Return true if this Data is lazy.
409      */
410      ESCRIPT_DLL_API
411      bool
412      isLazy() const;
413    
414      /**
415         \brief Return true if this data is ready.
416      */
417      ESCRIPT_DLL_API
418      bool
419      isReady() const;
420    
421      /**
422       \brief       \brief
423       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
424    contains datapoints.
425    */    */
426    ESCRIPT_DLL_API    ESCRIPT_DLL_API
427    bool    bool
# Line 395  class Data { Line 453  class Data {
453    */    */
454    ESCRIPT_DLL_API    ESCRIPT_DLL_API
455    inline    inline
456    const AbstractDomain&  //   const AbstractDomain&
457      const_Domain_ptr
458    getDomain() const    getDomain() const
459    {    {
460       return getFunctionSpace().getDomain();       return getFunctionSpace().getDomain();
461    }    }
462    
463    
464      /**
465         \brief
466         Return the domain.
467         TODO: For internal use only.   This should be removed.
468      */
469      ESCRIPT_DLL_API
470      inline
471    //   const AbstractDomain&
472      Domain_ptr
473      getDomainPython() const
474      {
475         return getFunctionSpace().getDomainPython();
476      }
477    
478    /**    /**
479       \brief       \brief
480       Return a copy of the domain.       Return a copy of the domain.
# Line 415  class Data { Line 489  class Data {
489    */    */
490    ESCRIPT_DLL_API    ESCRIPT_DLL_API
491    inline    inline
492    int    unsigned int
493    getDataPointRank() const    getDataPointRank() const
494    {    {
495      return m_data->getPointDataView().getRank();      return m_data->getRank();
496    }    }
497    
498    /**    /**
499       \brief       \brief
500         Return the number of data points
501      */
502      ESCRIPT_DLL_API
503      inline
504      int
505      getNumDataPoints() const
506      {
507        return getNumSamples() * getNumDataPointsPerSample();
508      }
509      /**
510         \brief
511       Return the number of samples.       Return the number of samples.
512    */    */
513    ESCRIPT_DLL_API    ESCRIPT_DLL_API
# Line 445  class Data { Line 530  class Data {
530      return m_data->getNumDPPSample();      return m_data->getNumDPPSample();
531    }    }
532    
533    
534      /**
535        \brief
536        Return the number of values in the shape for this object.
537      */
538      ESCRIPT_DLL_API
539      int
540      getNoValues() const
541      {
542        return m_data->getNoValues();
543      }
544    
545    
546      /**
547         \brief
548         dumps the object into a netCDF file
549      */
550      ESCRIPT_DLL_API
551      void
552      dump(const std::string fileName) const;
553    
554    /**    /**
555       \brief       \brief
556       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
# Line 454  class Data { Line 560  class Data {
560    ESCRIPT_DLL_API    ESCRIPT_DLL_API
561    inline    inline
562    DataAbstract::ValueType::value_type*    DataAbstract::ValueType::value_type*
563    getSampleData(DataAbstract::ValueType::size_type sampleNo)    getSampleData(DataAbstract::ValueType::size_type sampleNo);
   {  
     return m_data->getSampleData(sampleNo);  
   }  
564    
565    /**    /**
566       \brief       \brief
# Line 475  class Data { Line 578  class Data {
578    
579    /**    /**
580       \brief       \brief
581       Assign the given value to the data-points referenced by the given       Return a view into the data for the data point specified.
582       reference number.       NOTE: Construction of the DataArrayView is a relatively expensive
583         operation.
584       The value supplied is a python numarray object.  The data from this numarray       \param sampleNo - Input -
585       is unpacked into a DataArray, and this is used to set the corresponding       \param dataPointNo - Input -
      data-points in the underlying Data object.  
   
      If the underlying Data object cannot be accessed via reference numbers, an  
      exception will be thrown.  
   
      \param ref - Input - reference number.  
      \param value - Input - value to assign to data-points associated with  
                             the given reference number.  
586    */    */
587    ESCRIPT_DLL_API    ESCRIPT_DLL_API
588    void    DataTypes::ValueType::const_reference
589    setRefValue(int ref,    getDataPoint(int sampleNo, int dataPointNo) const;
               const boost::python::numeric::array& value);  
590    
   /**  
      \brief  
      Return the values associated with the data-points referenced by the given  
      reference number.  
591    
592       The value supplied is a python numarray object. The data from the corresponding    ESCRIPT_DLL_API
593       data-points in this Data object are packed into the given numarray object.    DataTypes::ValueType::reference
594      getDataPoint(int sampleNo, int dataPointNo);
595    
      If the underlying Data object cannot be accessed via reference numbers, an  
      exception will be thrown.  
596    
      \param ref - Input - reference number.  
      \param value - Output - object to receive values from data-points  
                              associated with the given reference number.  
   */  
   ESCRIPT_DLL_API  
   void  
   getRefValue(int ref,  
               boost::python::numeric::array& value);  
597    
598    /**    /**
599       \brief       \brief
600       Return a view into the data for the data point specified.       Return the offset for the given sample and point within the sample
      NOTE: Construction of the DataArrayView is a relatively expensive  
      operation.  
      \param sampleNo - Input -  
      \param dataPointNo - Input -  
601    */    */
602    ESCRIPT_DLL_API    ESCRIPT_DLL_API
603    inline    inline
604    DataArrayView    DataTypes::ValueType::size_type
605    getDataPoint(int sampleNo,    getDataOffset(int sampleNo,
606                 int dataPointNo)                 int dataPointNo)
607    {    {
608      return m_data->getDataPoint(sampleNo,dataPointNo);        return m_data->getPointOffset(sampleNo,dataPointNo);
609    }    }
610    
611    /**    /**
# Line 536  class Data { Line 613  class Data {
613       Return a reference to the data point shape.       Return a reference to the data point shape.
614    */    */
615    ESCRIPT_DLL_API    ESCRIPT_DLL_API
616    const DataArrayView::ShapeType&    inline
617    getDataPointShape() const;    const DataTypes::ShapeType&
618      getDataPointShape() const
619      {
620        return m_data->getShape();
621      }
622    
623    /**    /**
624       \brief       \brief
# Line 561  class Data { Line 642  class Data {
642       Return the number of doubles stored for this Data.       Return the number of doubles stored for this Data.
643    */    */
644    ESCRIPT_DLL_API    ESCRIPT_DLL_API
645    DataArrayView::ValueType::size_type    DataTypes::ValueType::size_type
646    getLength() const;    getLength() const;
647    
648    
649    
650    /**    /**
651       \brief       \brief
652       Assign the given value to the tag. Implicitly converts this       Assign the given value to the tag assocciated with name. Implicitly converts this
653       object to type DataTagged. Throws an exception if this object       object to type DataTagged. Throws an exception if this object
654       cannot be converted to a DataTagged object.       cannot be converted to a DataTagged object or name cannot be mapped onto a tag key.
655         \param tagKey - Input - Integer key.
656         \param value - Input - Value to associate with given key.
657        ==>*
658      */
659      ESCRIPT_DLL_API
660      void
661      setTaggedValueByName(std::string name,
662                           const boost::python::object& value);
663    
664      /**
665         \brief
666         Assign the given value to the tag. Implicitly converts this
667         object to type DataTagged if it is constant.
668    
669       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
670       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
671      ==>*      ==>*
# Line 581  class Data { Line 678  class Data {
678    /**    /**
679       \brief       \brief
680       Assign the given value to the tag. Implicitly converts this       Assign the given value to the tag. Implicitly converts this
681       object to type DataTagged. Throws an exception if this object       object to type DataTagged if it is constant.
682       cannot be converted to a DataTagged object.  
683       \param tagKey - Input - Integer key.       \param tagKey - Input - Integer key.
684         \param pointshape - Input - The shape of the value parameter
685       \param value - Input - Value to associate with given key.       \param value - Input - Value to associate with given key.
686      ==>*       \param dataOffset - Input - Offset of the begining of the point within the value parameter
687    */    */
688    ESCRIPT_DLL_API    ESCRIPT_DLL_API
689    void    void
690    setTaggedValueFromCPP(int tagKey,    setTaggedValueFromCPP(int tagKey,
691                          const DataArrayView& value);              const DataTypes::ShapeType& pointshape,
692                            const DataTypes::ValueType& value,
693                int dataOffset=0);
694    
695    
696    
697    /**    /**
698      \brief      \brief
# Line 607  class Data { Line 709  class Data {
709    
710    /**    /**
711       \brief       \brief
712         set all values to zero
713         *
714      */
715      ESCRIPT_DLL_API
716      void
717      setToZero();
718    
719      /**
720         \brief
721       Interpolates this onto the given functionspace and returns       Interpolates this onto the given functionspace and returns
722       the result as a Data object.       the result as a Data object.
723       *       *
# Line 614  class Data { Line 725  class Data {
725    ESCRIPT_DLL_API    ESCRIPT_DLL_API
726    Data    Data
727    interpolate(const FunctionSpace& functionspace) const;    interpolate(const FunctionSpace& functionspace) const;
   
728    /**    /**
729       \brief       \brief
730       Calculates the gradient of the data at the data points of functionspace.       Calculates the gradient of the data at the data points of functionspace.
# Line 636  class Data { Line 746  class Data {
746    */    */
747    ESCRIPT_DLL_API    ESCRIPT_DLL_API
748    boost::python::numeric::array    boost::python::numeric::array
749    integrate() const;    integrate_const() const;
750    
751      ESCRIPT_DLL_API
752      boost::python::numeric::array
753      integrate();
754    
755      /**
756         \brief
757         Returns 1./ Data object
758         *
759      */
760      ESCRIPT_DLL_API
761      Data
762      oneOver() const;
763    /**    /**
764       \brief       \brief
765       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.
# Line 695  class Data { Line 817  class Data {
817    /**    /**
818       \brief       \brief
819       Return the maximum absolute value of this Data object.       Return the maximum absolute value of this Data object.
820       *  
821         The method is not const because lazy data needs to be expanded before Lsup can be computed.
822         The _const form can be used when the Data object is const, however this will only work for
823         Data which is not Lazy.
824    
825         For Data which contain no samples (or tagged Data for which no tags in use have a value)
826         zero is returned.
827    */    */
828    ESCRIPT_DLL_API    ESCRIPT_DLL_API
829    double    double
830    Lsup() const;    Lsup();
831    
   /**  
      \brief  
      Return the minimum absolute value of this Data object.  
      *  
   */  
832    ESCRIPT_DLL_API    ESCRIPT_DLL_API
833    double    double
834    Linf() const;    Lsup_const() const;
835    
836    
837    /**    /**
838       \brief       \brief
839       Return the maximum value of this Data object.       Return the maximum value of this Data object.
840       *  
841         The method is not const because lazy data needs to be expanded before sup can be computed.
842         The _const form can be used when the Data object is const, however this will only work for
843         Data which is not Lazy.
844    
845         For Data which contain no samples (or tagged Data for which no tags in use have a value)
846         a large negative value is returned.
847    */    */
848    ESCRIPT_DLL_API    ESCRIPT_DLL_API
849    double    double
850    sup() const;    sup();
851    
852      ESCRIPT_DLL_API
853      double
854      sup_const() const;
855    
856    
857    /**    /**
858       \brief       \brief
859       Return the minimum value of this Data object.       Return the minimum value of this Data object.
860       *  
861         The method is not const because lazy data needs to be expanded before inf can be computed.
862         The _const form can be used when the Data object is const, however this will only work for
863         Data which is not Lazy.
864    
865         For Data which contain no samples (or tagged Data for which no tags in use have a value)
866         a large positive value is returned.
867    */    */
868    ESCRIPT_DLL_API    ESCRIPT_DLL_API
869    double    double
870    inf() const;    inf();
871    
872      ESCRIPT_DLL_API
873      double
874      inf_const() const;
875    
876    
877    
878    /**    /**
879       \brief       \brief
# Line 762  class Data { Line 909  class Data {
909    */    */
910    ESCRIPT_DLL_API    ESCRIPT_DLL_API
911    const boost::python::tuple    const boost::python::tuple
912    mindp() const;    minGlobalDataPoint() const;
913    
914    ESCRIPT_DLL_API    ESCRIPT_DLL_API
915    void    void
916    calc_mindp(int& SampleNo,    calc_minGlobalDataPoint(int& ProcNo,  int& DataPointNo) const;
              int& DataPointNo) const;  
   
917    /**    /**
918       \brief       \brief
919       Return the sign of each data point of this Data object.       Return the sign of each data point of this Data object.
# Line 781  class Data { Line 926  class Data {
926    
927    /**    /**
928       \brief       \brief
929         Return the symmetric part of a matrix which is half the matrix plus its transpose.
930         *
931      */
932      ESCRIPT_DLL_API
933      Data
934      symmetric() const;
935    
936      /**
937         \brief
938         Return the nonsymmetric part of a matrix which is half the matrix minus its transpose.
939         *
940      */
941      ESCRIPT_DLL_API
942      Data
943      nonsymmetric() const;
944    
945      /**
946         \brief
947         Return the trace of a matrix
948         *
949      */
950      ESCRIPT_DLL_API
951      Data
952      trace(int axis_offset) const;
953    
954      /**
955         \brief
956         Transpose each data point of this Data object around the given axis.
957         *
958      */
959      ESCRIPT_DLL_API
960      Data
961      transpose(int axis_offset) const;
962    
963      /**
964         \brief
965       Return the eigenvalues of the symmetric part at each data point of this Data object in increasing values.       Return the eigenvalues of the symmetric part at each data point of this Data object in increasing values.
966       Currently this function is restricted to rank 2, square shape, and dimension 3.       Currently this function is restricted to rank 2, square shape, and dimension 3.
967       *       *
# Line 792  class Data { Line 973  class Data {
973    /**    /**
974       \brief       \brief
975       Return the eigenvalues and corresponding eigenvcetors of the symmetric part at each data point of this Data object.       Return the eigenvalues and corresponding eigenvcetors of the symmetric part at each data point of this Data object.
976       the eigenvalues are ordered in increasing size where eigenvalues with relative difference less than       the eigenvalues are ordered in increasing size where eigenvalues with relative difference less than
977       tol are treated as equal. The eigenvectors are orthogonal, normalized and the sclaed such that the       tol are treated as equal. The eigenvectors are orthogonal, normalized and the sclaed such that the
978       first non-zero entry is positive.       first non-zero entry is positive.
979       Currently this function is restricted to rank 2, square shape, and dimension 3       Currently this function is restricted to rank 2, square shape, and dimension 3
980       *       *
# Line 804  class Data { Line 985  class Data {
985    
986    /**    /**
987       \brief       \brief
988       Transpose each data point of this Data object around the given axis.       swaps the components axis0 and axis1
      --* not implemented yet *--  
989       *       *
990    */    */
991    ESCRIPT_DLL_API    ESCRIPT_DLL_API
992    Data    Data
993    transpose(int axis) const;    swapaxes(const int axis0, const int axis1) const;
994    
995    /**    /**
996       \brief       \brief
997       Calculate the trace of each data point of this Data object.       Return the error function erf of each data point of this Data object.
998       *       *
999    */    */
1000    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1001    Data    Data
1002    trace() const;    erf() const;
1003    
1004    /**    /**
1005       \brief       \brief
# Line 998  class Data { Line 1178  class Data {
1178    /**    /**
1179       \brief       \brief
1180       Return the given power of each data point of this boost python object.       Return the given power of each data point of this boost python object.
1181        
1182       \param right Input - the power to raise the object to.       \param right Input - the power to raise the object to.
1183       *       *
1184     */     */
# Line 1009  class Data { Line 1189  class Data {
1189    /**    /**
1190       \brief       \brief
1191       Return the given power of each data point of this boost python object.       Return the given power of each data point of this boost python object.
1192        
1193       \param left Input - the bases       \param left Input - the bases
1194       *       *
1195     */     */
# Line 1045  class Data { Line 1225  class Data {
1225    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1226    Data& operator+=(const boost::python::object& right);    Data& operator+=(const boost::python::object& right);
1227    
1228      ESCRIPT_DLL_API
1229      Data& operator=(const Data& other);
1230    
1231    /**    /**
1232       \brief       \brief
1233       Overloaded operator -=       Overloaded operator -=
# Line 1137  class Data { Line 1320  class Data {
1320    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1321    inline    inline
1322    void    void
1323    unaryOp(UnaryFunction operation);    unaryOp2(UnaryFunction operation);
1324    
1325    /**    /**
1326       \brief       \brief
# Line 1148  class Data { Line 1331  class Data {
1331    */    */
1332    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1333    Data    Data
1334    getSlice(const DataArrayView::RegionType& region) const;    getSlice(const DataTypes::RegionType& region) const;
1335    
1336    /**    /**
1337       \brief       \brief
# Line 1161  class Data { Line 1344  class Data {
1344    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1345    void    void
1346    setSlice(const Data& value,    setSlice(const Data& value,
1347             const DataArrayView::RegionType& region);             const DataTypes::RegionType& region);
1348    
1349    /**    /**
1350       \brief       \brief
1351       Archive the current Data object to the given file.       print the data values to stdout. Used for debugging
      \param fileName - Input - file to archive to.  
1352    */    */
1353    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1354    void    void
1355    archiveData(const std::string fileName);          print(void);
1356    
1357    /**    /**
1358       \brief       \brief
1359       Extract the Data object archived in the given file, overwriting       return the MPI rank number of the local data
1360       the current Data object.                   MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()
1361       Note - the current object must be of type DataEmpty.                   is returned
      \param fileName - Input - file to extract from.  
      \param fspace - Input - a suitable FunctionSpace descibing the data.  
1362    */    */
1363    ESCRIPT_DLL_API    ESCRIPT_DLL_API
1364    void          int
1365    extractData(const std::string fileName,          get_MPIRank(void) const;
               const FunctionSpace& fspace);  
1366    
1367      /**
1368         \brief
1369         return the MPI rank number of the local data
1370                     MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()
1371                     is returned
1372      */
1373      ESCRIPT_DLL_API
1374            int
1375            get_MPISize(void) const;
1376    
1377    /**    /**
1378       \brief       \brief
1379       print the data values to stdout. Used for debugging       return the MPI rank number of the local data
1380                     MPI_COMM_WORLD is assumed and returned.
1381    */    */
1382    void print();    ESCRIPT_DLL_API
1383            MPI_Comm
1384            get_MPIComm(void) const;
1385    
1386      /**
1387         \brief
1388         return the object produced by the factory, which is a DataConstant or DataExpanded
1389        TODO Ownership of this object should be explained in doco.
1390      */
1391      ESCRIPT_DLL_API
1392            DataAbstract*
1393            borrowData(void) const;
1394    
1395      ESCRIPT_DLL_API
1396            DataAbstract_ptr
1397            borrowDataPtr(void) const;
1398    
1399      ESCRIPT_DLL_API
1400            DataReady_ptr
1401            borrowReadyPtr(void) const;
1402    
1403    
1404    
1405      /**
1406         \brief
1407         Return a pointer to the beginning of the datapoint at the specified offset.
1408         TODO Eventually these should be inlined.
1409         \param i - position(offset) in the underlying datastructure
1410      */
1411      ESCRIPT_DLL_API
1412            DataTypes::ValueType::const_reference
1413            getDataAtOffset(DataTypes::ValueType::size_type i) const;
1414    
1415    
1416      ESCRIPT_DLL_API
1417            DataTypes::ValueType::reference
1418            getDataAtOffset(DataTypes::ValueType::size_type i);
1419    
1420   protected:   protected:
1421    
1422   private:   private:
1423    
1424      double
1425      LsupWorker() const;
1426    
1427      double
1428      supWorker() const;
1429    
1430      double
1431      infWorker() const;
1432    
1433      boost::python::numeric::array
1434      integrateWorker() const;
1435    
1436    /**    /**
1437       \brief       \brief
1438       Check *this and the right operand are compatible. Throws       Check *this and the right operand are compatible. Throws
# Line 1249  class Data { Line 1486  class Data {
1486    
1487    /**    /**
1488       \brief       \brief
      Perform the given binary operation on all of the data's elements.  
      RHS is a boost::python object.  
   */  
   template <class BinaryFunction>  
   inline  
   void  
   binaryOp(const boost::python::object& right,  
            BinaryFunction operation);  
   
   /**  
      \brief  
1489       Convert the data type of the RHS to match this.       Convert the data type of the RHS to match this.
1490       \param right - Input - data type to match.       \param right - Input - data type to match.
1491    */    */
# Line 1278  class Data { Line 1504  class Data {
1504       \brief       \brief
1505       Construct a Data object of the appropriate type.       Construct a Data object of the appropriate type.
1506    */    */
1507    template <class IValueType>  
1508    void    void
1509    initialise(const IValueType& value,    initialise(const DataTypes::ValueType& value,
1510             const DataTypes::ShapeType& shape,
1511               const FunctionSpace& what,               const FunctionSpace& what,
1512               bool expanded);               bool expanded);
1513    
   /**  
      \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.  
   */  
1514    void    void
1515    reshapeDataPoint(const DataArrayView::ShapeType& shape);    initialise(const boost::python::numeric::array& value,
1516                     const FunctionSpace& what,
1517                     bool expanded);
1518    
1519    //    //
1520    // pointer to the actual data object    // flag to protect the data object against any update
1521    boost::shared_ptr<DataAbstract> m_data;    bool m_protected;
1522    
1523    //    //
1524    // pointer to the internal profiling data    // pointer to the actual data object
1525    struct profDataEntry *profData;  //   boost::shared_ptr<DataAbstract> m_data;
1526      DataAbstract_ptr m_data;
1527    
1528    // If possible please use getReadyPtr instead
1529      const DataReady*
1530      getReady() const;
1531    
1532      DataReady*
1533      getReady();
1534    
1535      DataReady_ptr
1536      getReadyPtr();
1537    
1538      const_DataReady_ptr
1539      getReadyPtr() const;
1540    
1541    
1542  };  };
1543    
1544  template <class IValueType>  }   // end namespace escript
1545  void  
1546  Data::initialise(const IValueType& value,  
1547                   const FunctionSpace& what,  // No, this is not supposed to be at the top of the file
1548                   bool expanded)  // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1549    // so that I can dynamic cast between them below.
1550    #include "DataReady.h"
1551    
1552    namespace escript
1553  {  {
1554    //  
1555    // Construct a Data object of the appropriate type.  inline
1556    // Construct the object first as there seems to be a bug which causes  const DataReady*
1557    // undefined behaviour if an exception is thrown during construction  Data::getReady() const
1558    // within the shared_ptr constructor.  {
1559    if (expanded) {     const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1560      DataAbstract* temp=new DataExpanded(value,what);     EsysAssert((dr!=0), "Error - casting to DataReady.");
1561      boost::shared_ptr<DataAbstract> temp_data(temp);     return dr;
1562      m_data=temp_data;  }
1563    } else {  
1564      DataAbstract* temp=new DataConstant(value,what);  inline
1565      boost::shared_ptr<DataAbstract> temp_data(temp);  DataReady*
1566      m_data=temp_data;  Data::getReady()
1567    }  {
1568       DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1569       EsysAssert((dr!=0), "Error - casting to DataReady.");
1570       return dr;
1571    }
1572    
1573    inline
1574    DataReady_ptr
1575    Data::getReadyPtr()
1576    {
1577       DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1578       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1579       return dr;
1580  }  }
1581    
1582    
1583    inline
1584    const_DataReady_ptr
1585    Data::getReadyPtr() const
1586    {
1587       const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1588       EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1589       return dr;
1590    }
1591    
1592    inline
1593    DataAbstract::ValueType::value_type*
1594    Data::getSampleData(DataAbstract::ValueType::size_type sampleNo)
1595    {
1596       if (isLazy())
1597       {
1598        resolve();
1599       }
1600       return getReady()->getSampleData(sampleNo);
1601    }
1602    
1603    
1604    /**
1605       Modify a filename for MPI parallel output to multiple files
1606    */
1607    char *Escript_MPI_appendRankToFileName(const char *, int, int);
1608    
1609  /**  /**
1610     Binary Data object operators.     Binary Data object operators.
1611  */  */
1612    inline double rpow(double x,double y)
1613    {
1614        return pow(y,x);
1615    }
1616    
1617  /**  /**
1618    \brief    \brief
# Line 1422  ESCRIPT_DLL_API Data operator*(const boo Line 1706  ESCRIPT_DLL_API Data operator*(const boo
1706  */  */
1707  ESCRIPT_DLL_API Data operator/(const boost::python::object& left, const Data& right);  ESCRIPT_DLL_API Data operator/(const boost::python::object& left, const Data& right);
1708    
1709    
1710    
1711  /**  /**
1712    \brief    \brief
1713    Output operator    Output operator
# Line 1430  ESCRIPT_DLL_API std::ostream& operator<< Line 1716  ESCRIPT_DLL_API std::ostream& operator<<
1716    
1717  /**  /**
1718    \brief    \brief
1719    Return true if operands are equivalent, else return false.    Compute a tensor product of two Data objects
1720    NB: this operator does very little at this point, and isn't to    \param arg0 - Input - Data object
1721    be relied on. Requires further implementation.    \param arg1 - Input - Data object
1722      \param axis_offset - Input - axis offset
1723      \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1724  */  */
1725  //ESCRIPT_DLL_API bool operator==(const Data& left, const Data& right);  ESCRIPT_DLL_API
1726    Data
1727    C_GeneralTensorProduct(Data& arg0,
1728                         Data& arg1,
1729                         int axis_offset=0,
1730                         int transpose=0);
1731    
1732  /**  /**
1733    \brief    \brief
# Line 1449  Data::binaryOp(const Data& right, Line 1742  Data::binaryOp(const Data& right,
1742  {  {
1743     //     //
1744     // 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
1745     if (getPointDataView().getRank()==0 && right.getPointDataView().getRank()!=0) {     if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1746       reshapeDataPoint(right.getPointDataView().getShape());       throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
1747       }
1748    
1749       if (isLazy() || right.isLazy())
1750       {
1751         throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
1752     }     }
1753     //     //
1754     // initially make the temporary a shallow copy     // initially make the temporary a shallow copy
1755     Data tempRight(right);     Data tempRight(right);
1756    
1757     if (getFunctionSpace()!=right.getFunctionSpace()) {     if (getFunctionSpace()!=right.getFunctionSpace()) {
1758       if (right.probeInterpolation(getFunctionSpace())) {       if (right.probeInterpolation(getFunctionSpace())) {
1759         //         //
1760         // an interpolation is required so create a new Data         // an interpolation is required so create a new Data
1761         tempRight=Data(right,this->getFunctionSpace());         tempRight=Data(right,this->getFunctionSpace());
1762       } else if (probeInterpolation(right.getFunctionSpace())) {       } else if (probeInterpolation(right.getFunctionSpace())) {
1763         //         //
# Line 1480  Data::binaryOp(const Data& right, Line 1779  Data::binaryOp(const Data& right,
1779       // of any data type       // of any data type
1780       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());       DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
1781       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");       EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
1782       escript::binaryOp(*leftC,*(tempRight.m_data.get()),operation);       escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
1783     } else if (isTagged()) {     } else if (isTagged()) {
1784       //       //
1785       // 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 1506  Data::binaryOp(const Data& right, Line 1805  Data::binaryOp(const Data& right,
1805    
1806  /**  /**
1807    \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  
1808    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.
1809    Given operation combines each element of each data point, thus argument    Given operation combines each element of each data point, thus argument
1810    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.
# Line 1614  Data::algorithm(BinaryFunction operation Line 1827  Data::algorithm(BinaryFunction operation
1827      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
1828      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");      EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
1829      return escript::algorithm(*leftC,operation,initial_value);      return escript::algorithm(*leftC,operation,initial_value);
1830      } else if (isEmpty()) {
1831        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1832      } else if (isLazy()) {
1833        throw DataException("Error - Operations not permitted on instances of DataLazy.");
1834      } else {
1835        throw DataException("Error - Data encapsulates an unknown type.");
1836    }    }
   return 0;  
1837  }  }
1838    
1839  /**  /**
1840    \brief    \brief
1841    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.
1842    Given operation combines each element within each data point into a scalar,    Given operation combines each element within each data point into a scalar,
1843    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
1844    rank 0 Data object.    rank 0 Data object.
1845    Calls escript::dp_algorithm.    Calls escript::dp_algorithm.
1846  */  */
# Line 1631  inline Line 1849  inline
1849  Data  Data
1850  Data::dp_algorithm(BinaryFunction operation, double initial_value) const  Data::dp_algorithm(BinaryFunction operation, double initial_value) const
1851  {  {
1852    if (isExpanded()) {    if (isEmpty()) {
1853      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());      throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1854      }
1855      else if (isExpanded()) {
1856        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
1857      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
1858      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());      DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
1859      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");      EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
1860      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");      EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
1861      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);      escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
1862      return result;      return result;
1863    } else if (isTagged()) {    }
1864      else if (isTagged()) {
1865      DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());
     DataArrayView::ShapeType viewShape;  
     DataArrayView::ValueType viewData(1);  
     viewData[0]=0;  
     DataArrayView defaultValue(viewData,viewShape);  
     DataTagged::TagListType keys;  
     DataTagged::ValueListType values;  
     DataTagged::DataMapType::const_iterator i;  
     for (i=dataT->getTagLookup().begin();i!=dataT->getTagLookup().end();i++) {  
       keys.push_back(i->first);  
       values.push_back(defaultValue);  
     }  
     Data result(keys,values,defaultValue,getFunctionSpace());  
     DataTagged* resultT=dynamic_cast<DataTagged*>(result.m_data.get());  
1866      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");      EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
1867      EsysAssert((resultT!=0), "Programming error - casting result to DataTagged.");      DataTypes::ValueType defval(1);
1868        defval[0]=0;
1869        DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
1870      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);      escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
1871      return result;      return Data(resultT);   // note: the Data object now owns the resultT pointer
1872    } else if (isConstant()) {    }
1873      Data result(0,DataArrayView::ShapeType(),getFunctionSpace(),isExpanded());    else if (isConstant()) {
1874        Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
1875      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
1876      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());      DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
1877      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");      EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
1878      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");      EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
1879      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);      escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
1880      return result;      return result;
1881      } else if (isLazy()) {
1882        throw DataException("Error - Operations not permitted on instances of DataLazy.");
1883      } else {
1884        throw DataException("Error - Data encapsulates an unknown type.");
1885    }    }
1886    Data falseRetVal; // to keep compiler quiet  }
1887    return falseRetVal;  
1888    /**
1889      \brief
1890      Compute a tensor operation with two Data objects
1891      \param arg0 - Input - Data object
1892      \param arg1 - Input - Data object
1893      \param operation - Input - Binary op functor
1894    */
1895    template <typename BinaryFunction>
1896    inline
1897    Data
1898    C_TensorBinaryOperation(Data const &arg_0,
1899                            Data const &arg_1,
1900                            BinaryFunction operation)
1901    {
1902      if (arg_0.isEmpty() || arg_1.isEmpty())
1903      {
1904         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1905      }
1906      if (arg_0.isLazy() || arg_1.isLazy())
1907      {
1908         throw DataException("Error - Operations not permitted on lazy data.");
1909      }
1910      // Interpolate if necessary and find an appropriate function space
1911      Data arg_0_Z, arg_1_Z;
1912      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
1913        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
1914          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
1915          arg_1_Z = Data(arg_1);
1916        }
1917        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
1918          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
1919          arg_0_Z =Data(arg_0);
1920        }
1921        else {
1922          throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible function spaces.");
1923        }
1924      } else {
1925          arg_0_Z = Data(arg_0);
1926          arg_1_Z = Data(arg_1);
1927      }
1928      // Get rank and shape of inputs
1929      int rank0 = arg_0_Z.getDataPointRank();
1930      int rank1 = arg_1_Z.getDataPointRank();
1931      DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
1932      DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
1933      int size0 = arg_0_Z.getDataPointSize();
1934      int size1 = arg_1_Z.getDataPointSize();
1935    
1936      // Declare output Data object
1937      Data res;
1938    
1939      if (shape0 == shape1) {
1940    
1941        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
1942          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
1943    /*      double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
1944          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
1945          double *ptr_2 = &((res.getPointDataView().getData())[0]);*/
1946          double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));
1947          double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));
1948          double *ptr_2 = &(res.getDataAtOffset(0));
1949    
1950          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
1951        }
1952        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
1953    
1954          // Prepare the DataConstant input
1955          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
1956    
1957          // Borrow DataTagged input from Data object
1958          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
1959    
1960          // Prepare a DataTagged output 2
1961          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
1962          res.tag();
1963          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
1964    
1965          // Prepare offset into DataConstant
1966          int offset_0 = tmp_0->getPointOffset(0,0);
1967          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
1968          // Get the views
1969    //       DataArrayView view_1 = tmp_1->getDefaultValue();
1970    //       DataArrayView view_2 = tmp_2->getDefaultValue();
1971    //       // Get the pointers to the actual data
1972    //       double *ptr_1 = &((view_1.getData())[0]);
1973    //       double *ptr_2 = &((view_2.getData())[0]);
1974    
1975          // Get the pointers to the actual data
1976          double *ptr_1 = &(tmp_1->getDefaultValue(0));
1977          double *ptr_2 = &(tmp_2->getDefaultValue(0));
1978    
1979          // Compute a result for the default
1980          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
1981          // Compute a result for each tag
1982          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
1983          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
1984          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
1985            tmp_2->addTag(i->first);
1986    /*        DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
1987            DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
1988            double *ptr_1 = &view_1.getData(0);
1989            double *ptr_2 = &view_2.getData(0);*/
1990            double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
1991            double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
1992    
1993            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
1994          }
1995    
1996        }
1997        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
1998    
1999          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2000          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2001          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2002          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2003    
2004          int sampleNo_1,dataPointNo_1;
2005          int numSamples_1 = arg_1_Z.getNumSamples();
2006          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2007          int offset_0 = tmp_0->getPointOffset(0,0);
2008          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2009          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2010            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2011              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2012              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2013    //           double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2014    //           double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2015    //           double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2016    
2017              double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2018              double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2019              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2020              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2021            }
2022          }
2023    
2024        }
2025        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2026    
2027          // Borrow DataTagged input from Data object
2028          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2029    
2030          // Prepare the DataConstant input
2031          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2032    
2033          // Prepare a DataTagged output 2
2034          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2035          res.tag();
2036          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2037    
2038          // Prepare offset into DataConstant
2039          int offset_1 = tmp_1->getPointOffset(0,0);
2040    //       double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2041          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2042          // Get the views
2043    //       DataArrayView view_0 = tmp_0->getDefaultValue();
2044    //       DataArrayView view_2 = tmp_2->getDefaultValue();
2045    //       // Get the pointers to the actual data
2046    //       double *ptr_0 = &((view_0.getData())[0]);
2047    //       double *ptr_2 = &((view_2.getData())[0]);
2048          // Get the pointers to the actual data
2049          double *ptr_0 = &(tmp_0->getDefaultValue(0));
2050          double *ptr_2 = &(tmp_2->getDefaultValue(0));
2051          // Compute a result for the default
2052          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2053          // Compute a result for each tag
2054          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2055          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2056          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2057            tmp_2->addTag(i->first);
2058    //         DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2059    //         DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2060    //         double *ptr_0 = &view_0.getData(0);
2061    //         double *ptr_2 = &view_2.getData(0);
2062            double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2063            double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2064            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2065          }
2066    
2067        }
2068        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2069    
2070          // Borrow DataTagged input from Data object
2071          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2072    
2073          // Borrow DataTagged input from Data object
2074          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2075    
2076          // Prepare a DataTagged output 2
2077          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2078          res.tag();        // DataTagged output
2079          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2080    
2081    //       // Get the views
2082    //       DataArrayView view_0 = tmp_0->getDefaultValue();
2083    //       DataArrayView view_1 = tmp_1->getDefaultValue();
2084    //       DataArrayView view_2 = tmp_2->getDefaultValue();
2085    //       // Get the pointers to the actual data
2086    //       double *ptr_0 = &((view_0.getData())[0]);
2087    //       double *ptr_1 = &((view_1.getData())[0]);
2088    //       double *ptr_2 = &((view_2.getData())[0]);
2089    
2090          // Get the pointers to the actual data
2091          double *ptr_0 = &(tmp_0->getDefaultValue(0));
2092          double *ptr_1 = &(tmp_1->getDefaultValue(0));
2093          double *ptr_2 = &(tmp_2->getDefaultValue(0));
2094    
2095          // Compute a result for the default
2096          tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2097          // Merge the tags
2098          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2099          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2100          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2101          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2102            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2103          }
2104          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2105            tmp_2->addTag(i->first);
2106          }
2107          // Compute a result for each tag
2108          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2109          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2110    
2111    //         DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2112    //         DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2113    //         DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2114    //         double *ptr_0 = &view_0.getData(0);
2115    //         double *ptr_1 = &view_1.getData(0);
2116    //         double *ptr_2 = &view_2.getData(0);
2117    
2118            double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2119            double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2120            double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2121    
2122            tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2123          }
2124    
2125        }
2126        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2127    
2128          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2129          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2130          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2131          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2132          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2133    
2134          int sampleNo_0,dataPointNo_0;
2135          int numSamples_0 = arg_0_Z.getNumSamples();
2136          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2137          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2138          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2139            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2140            double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2141            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2142              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2143              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2144    
2145    //           double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2146    //           double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2147              double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2148              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2149    
2150    
2151              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2152            }
2153          }
2154    
2155        }
2156        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2157    
2158          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2159          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2160          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2161          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2162    
2163          int sampleNo_0,dataPointNo_0;
2164          int numSamples_0 = arg_0_Z.getNumSamples();
2165          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2166          int offset_1 = tmp_1->getPointOffset(0,0);
2167          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2168          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2169            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2170              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2171              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2172    
2173    //           double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2174    //           double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2175    //           double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2176    
2177              double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2178              double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2179              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2180    
2181    
2182              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2183            }
2184          }
2185    
2186        }
2187        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2188    
2189          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2190          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2191          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2192          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2193          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2194    
2195          int sampleNo_0,dataPointNo_0;
2196          int numSamples_0 = arg_0_Z.getNumSamples();
2197          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2198          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2199          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2200            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2201            double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2202            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2203              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2204              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2205              double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2206              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2207              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2208            }
2209          }
2210    
2211        }
2212        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2213    
2214          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2215          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2216          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2217          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2218          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2219    
2220          int sampleNo_0,dataPointNo_0;
2221          int numSamples_0 = arg_0_Z.getNumSamples();
2222          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2223          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2224          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2225            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2226              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2227              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2228              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2229              double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2230              double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2231              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2232              tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2233            }
2234          }
2235    
2236        }
2237        else {
2238          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2239        }
2240    
2241      } else if (0 == rank0) {
2242    
2243        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2244          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataConstant output
2245          double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));
2246          double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));
2247          double *ptr_2 = &(res.getDataAtOffset(0));
2248          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2249        }
2250        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2251    
2252          // Prepare the DataConstant input
2253          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2254    
2255          // Borrow DataTagged input from Data object
2256          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2257    
2258          // Prepare a DataTagged output 2
2259          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());      // DataTagged output
2260          res.tag();
2261          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2262    
2263          // Prepare offset into DataConstant
2264          int offset_0 = tmp_0->getPointOffset(0,0);
2265          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2266          // Get the views
2267    //       DataArrayView view_1 = tmp_1->getDefaultValue();
2268    //       DataArrayView view_2 = tmp_2->getDefaultValue();
2269    //       // Get the pointers to the actual data
2270    //       double *ptr_1 = &((view_1.getData())[0]);
2271    //       double *ptr_2 = &((view_2.getData())[0]);
2272           double *ptr_1 = &(tmp_1->getDefaultValue(0));
2273           double *ptr_2 = &(tmp_2->getDefaultValue(0));
2274    
2275          // Compute a result for the default
2276          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2277          // Compute a result for each tag
2278          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2279          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2280          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2281            tmp_2->addTag(i->first);
2282    //         DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2283    //         DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2284    //         double *ptr_1 = &view_1.getData(0);
2285    //         double *ptr_2 = &view_2.getData(0);
2286            double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2287            double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2288            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2289          }
2290    
2291        }
2292        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2293    
2294          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2295          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2296          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2297          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2298    
2299          int sampleNo_1,dataPointNo_1;
2300          int numSamples_1 = arg_1_Z.getNumSamples();
2301          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2302          int offset_0 = tmp_0->getPointOffset(0,0);
2303          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2304          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2305            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2306              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2307              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2308              double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2309              double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2310              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2311              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2312    
2313            }
2314          }
2315    
2316        }
2317        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2318    
2319          // Borrow DataTagged input from Data object
2320          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2321    
2322          // Prepare the DataConstant input
2323          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2324    
2325          // Prepare a DataTagged output 2
2326          res = Data(0.0, shape1, arg_0_Z.getFunctionSpace());      // DataTagged output
2327          res.tag();
2328          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2329    
2330          // Prepare offset into DataConstant
2331          int offset_1 = tmp_1->getPointOffset(0,0);
2332    //       double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2333          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2334          // Get the views
2335    /*      DataArrayView view_0 = tmp_0->getDefaultValue();
2336          DataArrayView view_2 = tmp_2->getDefaultValue();
2337          // Get the pointers to the actual data
2338          double *ptr_0 = &((view_0.getData())[0]);
2339          double *ptr_2 = &((view_2.getData())[0]);*/
2340    
2341          // Get the pointers to the actual data
2342          double *ptr_0 = &(tmp_0->getDefaultValue(0));
2343          double *ptr_2 = &(tmp_2->getDefaultValue(0));
2344    
2345    
2346          // Compute a result for the default
2347          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2348          // Compute a result for each tag
2349          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2350          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2351          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2352            tmp_2->addTag(i->first);
2353    /*        DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2354            DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2355            double *ptr_0 = &view_0.getData(0);
2356            double *ptr_2 = &view_2.getData(0);*/
2357            double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2358            double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2359    
2360            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2361          }
2362    
2363        }
2364        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2365    
2366          // Borrow DataTagged input from Data object
2367          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2368    
2369          // Borrow DataTagged input from Data object
2370          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2371    
2372          // Prepare a DataTagged output 2
2373          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
2374          res.tag();        // DataTagged output
2375          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2376    
2377          // Get the views
2378    /*      DataArrayView view_0 = tmp_0->getDefaultValue();
2379          DataArrayView view_1 = tmp_1->getDefaultValue();
2380          DataArrayView view_2 = tmp_2->getDefaultValue();
2381          // Get the pointers to the actual data
2382          double *ptr_0 = &((view_0.getData())[0]);
2383          double *ptr_1 = &((view_1.getData())[0]);
2384          double *ptr_2 = &((view_2.getData())[0]);*/
2385    
2386          // Get the pointers to the actual data
2387          double *ptr_0 = &(tmp_0->getDefaultValue(0));
2388          double *ptr_1 = &(tmp_1->getDefaultValue(0));
2389          double *ptr_2 = &(tmp_2->getDefaultValue(0));
2390    
2391    
2392          // Compute a result for the default
2393          tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2394          // Merge the tags
2395          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2396          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2397          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2398          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2399            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2400          }
2401          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2402            tmp_2->addTag(i->first);
2403          }
2404          // Compute a result for each tag
2405          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2406          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2407    
2408    /*        DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2409            DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2410            DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2411            double *ptr_0 = &view_0.getData(0);
2412            double *ptr_1 = &view_1.getData(0);
2413            double *ptr_2 = &view_2.getData(0);*/
2414    
2415            double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2416            double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2417            double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2418    
2419            tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2420          }
2421    
2422        }
2423        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2424    
2425          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2426          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2427          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2428          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2429          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2430    
2431          int sampleNo_0,dataPointNo_0;
2432          int numSamples_0 = arg_0_Z.getNumSamples();
2433          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2434          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2435          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2436            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2437            double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2438            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2439              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2440              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2441              double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2442              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2443              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2444            }
2445          }
2446    
2447        }
2448        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2449    
2450          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2451          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2452          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2453          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2454    
2455          int sampleNo_0,dataPointNo_0;
2456          int numSamples_0 = arg_0_Z.getNumSamples();
2457          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2458          int offset_1 = tmp_1->getPointOffset(0,0);
2459          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2460          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2461            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2462              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2463              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2464              double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2465              double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2466              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2467              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2468            }
2469          }
2470    
2471    
2472        }
2473        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2474    
2475          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2476          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2477          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2478          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2479          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2480    
2481          int sampleNo_0,dataPointNo_0;
2482          int numSamples_0 = arg_0_Z.getNumSamples();
2483          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2484          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2485          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2486            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2487            double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2488            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2489              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2490              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2491              double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2492              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2493              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2494            }
2495          }
2496    
2497        }
2498        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2499    
2500          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2501          res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2502          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2503          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2504          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2505    
2506          int sampleNo_0,dataPointNo_0;
2507          int numSamples_0 = arg_0_Z.getNumSamples();
2508          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2509          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2510          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2511            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2512              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2513              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2514              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2515              double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2516              double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2517              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2518              tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2519            }
2520          }
2521    
2522        }
2523        else {
2524          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2525        }
2526    
2527      } else if (0 == rank1) {
2528    
2529        if (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2530          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataConstant output
2531          double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));
2532          double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));
2533          double *ptr_2 = &(res.getDataAtOffset(0));
2534          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2535        }
2536        else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2537    
2538          // Prepare the DataConstant input
2539          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2540    
2541          // Borrow DataTagged input from Data object
2542          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2543    
2544          // Prepare a DataTagged output 2
2545          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());      // DataTagged output
2546          res.tag();
2547          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2548    
2549          // Prepare offset into DataConstant
2550          int offset_0 = tmp_0->getPointOffset(0,0);
2551          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2552          // Get the views
2553    /*      DataArrayView view_1 = tmp_1->getDefaultValue();
2554          DataArrayView view_2 = tmp_2->getDefaultValue();
2555          // Get the pointers to the actual data
2556          double *ptr_1 = &((view_1.getData())[0]);
2557          double *ptr_2 = &((view_2.getData())[0]);*/
2558          //Get the pointers to the actual data
2559          double *ptr_1 = &(tmp_1->getDefaultValue(0));
2560          double *ptr_2 = &(tmp_2->getDefaultValue(0));
2561    
2562          // Compute a result for the default
2563          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2564          // Compute a result for each tag
2565          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2566          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2567          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2568            tmp_2->addTag(i->first);
2569    //         DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2570    //         DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2571    //         double *ptr_1 = &view_1.getData(0);
2572    //         double *ptr_2 = &view_2.getData(0);
2573            double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2574            double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2575    
2576    
2577            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2578          }
2579    
2580        }
2581        else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2582    
2583          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2584          DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2585          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2586          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2587    
2588          int sampleNo_1,dataPointNo_1;
2589          int numSamples_1 = arg_1_Z.getNumSamples();
2590          int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2591          int offset_0 = tmp_0->getPointOffset(0,0);
2592          #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2593          for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2594            for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2595              int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2596              int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2597              double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2598              double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2599              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2600              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2601            }
2602          }
2603    
2604        }
2605        else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2606    
2607          // Borrow DataTagged input from Data object
2608          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2609    
2610          // Prepare the DataConstant input
2611          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2612    
2613          // Prepare a DataTagged output 2
2614          res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataTagged output
2615          res.tag();
2616          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2617    
2618          // Prepare offset into DataConstant
2619          int offset_1 = tmp_1->getPointOffset(0,0);
2620          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2621          // Get the views
2622    //       DataArrayView view_0 = tmp_0->getDefaultValue();
2623    //       DataArrayView view_2 = tmp_2->getDefaultValue();
2624    //       // Get the pointers to the actual data
2625    //       double *ptr_0 = &((view_0.getData())[0]);
2626    //       double *ptr_2 = &((view_2.getData())[0]);
2627          // Get the pointers to the actual data
2628          double *ptr_0 = &(tmp_0->getDefaultValue(0));
2629          double *ptr_2 = &(tmp_2->getDefaultValue(0));
2630          // Compute a result for the default
2631          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2632          // Compute a result for each tag
2633          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2634          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2635          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2636            tmp_2->addTag(i->first);
2637    /*        DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2638            DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2639            double *ptr_0 = &view_0.getData(0);
2640            double *ptr_2 = &view_2.getData(0);*/
2641            double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2642            double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2643            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2644          }
2645    
2646        }
2647        else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2648    
2649          // Borrow DataTagged input from Data object
2650          DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2651    
2652          // Borrow DataTagged input from Data object
2653          DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2654    
2655          // Prepare a DataTagged output 2
2656          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2657          res.tag();        // DataTagged output
2658          DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2659    
2660          // Get the views
2661    //       DataArrayView view_0 = tmp_0->getDefaultValue();
2662    //       DataArrayView view_1 = tmp_1->getDefaultValue();
2663    //       DataArrayView view_2 = tmp_2->getDefaultValue();
2664    //       // Get the pointers to the actual data
2665    //       double *ptr_0 = &((view_0.getData())[0]);
2666    //       double *ptr_1 = &((view_1.getData())[0]);
2667    //       double *ptr_2 = &((view_2.getData())[0]);
2668    
2669          // Get the pointers to the actual data
2670          double *ptr_0 = &(tmp_0->getDefaultValue(0));
2671          double *ptr_1 = &(tmp_1->getDefaultValue(0));
2672          double *ptr_2 = &(tmp_2->getDefaultValue(0));
2673    
2674          // Compute a result for the default
2675          tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2676          // Merge the tags
2677          DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2678          const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2679          const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2680          for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2681            tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2682          }
2683          for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2684            tmp_2->addTag(i->first);
2685          }
2686          // Compute a result for each tag
2687          const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2688          for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2689    //         DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2690    //         DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2691    //         DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2692    //         double *ptr_0 = &view_0.getData(0);
2693    //         double *ptr_1 = &view_1.getData(0);
2694    //         double *ptr_2 = &view_2.getData(0);
2695    
2696            double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2697            double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2698            double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2699            tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2700          }
2701    
2702        }
2703        else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2704    
2705          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2706          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2707          DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2708          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2709          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2710    
2711          int sampleNo_0,dataPointNo_0;
2712          int numSamples_0 = arg_0_Z.getNumSamples();
2713          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2714          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2715          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2716            int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2717            double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2718            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2719              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2720              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2721              double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2722              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2723              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2724            }
2725          }
2726    
2727        }
2728        else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2729    
2730          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2731          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2732          DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2733          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2734    
2735          int sampleNo_0,dataPointNo_0;
2736          int numSamples_0 = arg_0_Z.getNumSamples();
2737          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2738          int offset_1 = tmp_1->getPointOffset(0,0);
2739          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2740          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2741            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2742              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2743              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2744              double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2745              double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2746              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2747              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2748            }
2749          }
2750    
2751    
2752        }
2753        else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2754    
2755          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2756          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2757          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2758          DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2759          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2760    
2761          int sampleNo_0,dataPointNo_0;
2762          int numSamples_0 = arg_0_Z.getNumSamples();
2763          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2764          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2765          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2766            int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2767            double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2768            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2769              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2770              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2771              double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2772              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2773              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2774            }
2775          }
2776    
2777        }
2778        else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2779    
2780          // After finding a common function space above the two inputs have the same numSamples and num DPPS
2781          res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2782          DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2783          DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2784          DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2785    
2786          int sampleNo_0,dataPointNo_0;
2787          int numSamples_0 = arg_0_Z.getNumSamples();
2788          int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2789          #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2790          for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2791            for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2792              int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2793              int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2794              int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2795              double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2796              double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2797              double *ptr_2 = &(res.getDataAtOffset(offset_2));
2798              tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2799            }
2800          }
2801    
2802        }
2803        else {
2804          throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2805        }
2806    
2807      } else {
2808        throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible shapes");
2809      }
2810    
2811      return res;
2812    }
2813    
2814    template <typename UnaryFunction>
2815    Data
2816    C_TensorUnaryOperation(Data const &arg_0,
2817                           UnaryFunction operation)
2818    {
2819      if (arg_0.isEmpty())  // do this before we attempt to interpolate
2820      {
2821         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2822      }
2823      if (arg_0.isLazy())
2824      {
2825         throw DataException("Error - Operations not permitted on lazy data.");
2826      }
2827      // Interpolate if necessary and find an appropriate function space
2828      Data arg_0_Z = Data(arg_0);
2829    
2830      // Get rank and shape of inputs
2831    //  int rank0 = arg_0_Z.getDataPointRank();
2832      const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2833      int size0 = arg_0_Z.getDataPointSize();
2834    
2835      // Declare output Data object
2836      Data res;
2837    
2838      if (arg_0_Z.isConstant()) {
2839        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());      // DataConstant output
2840    //     double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2841    //     double *ptr_2 = &((res.getPointDataView().getData())[0]);
2842        double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));
2843        double *ptr_2 = &(res.getDataAtOffset(0));
2844        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2845      }
2846      else if (arg_0_Z.isTagged()) {
2847    
2848        // Borrow DataTagged input from Data object
2849        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2850    
2851        // Prepare a DataTagged output 2
2852        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace());   // DataTagged output
2853        res.tag();
2854        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2855    
2856    //     // Get the views
2857    //     DataArrayView view_0 = tmp_0->getDefaultValue();
2858    //     DataArrayView view_2 = tmp_2->getDefaultValue();
2859    //     // Get the pointers to the actual data
2860    //     double *ptr_0 = &((view_0.getData())[0]);
2861    //     double *ptr_2 = &((view_2.getData())[0]);
2862        // Get the pointers to the actual data
2863        double *ptr_0 = &(tmp_0->getDefaultValue(0));
2864        double *ptr_2 = &(tmp_2->getDefaultValue(0));
2865        // Compute a result for the default
2866        tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2867        // Compute a result for each tag
2868        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2869        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2870        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2871          tmp_2->addTag(i->first);
2872    //       DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2873    //       DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2874    //       double *ptr_0 = &view_0.getData(0);
2875    //       double *ptr_2 = &view_2.getData(0);
2876          double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2877          double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2878          tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2879        }
2880    
2881      }
2882      else if (arg_0_Z.isExpanded()) {
2883    
2884        res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
2885        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2886        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2887    
2888        int sampleNo_0,dataPointNo_0;
2889        int numSamples_0 = arg_0_Z.getNumSamples();
2890        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2891        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2892        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2893          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2894    //         int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2895    //         int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2896    //         double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2897    //         double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2898            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2899            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2900            double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2901            double *ptr_2 = &(res.getDataAtOffset(offset_2));
2902            tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2903          }
2904        }
2905      }
2906      else {
2907        throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
2908      }
2909    
2910      return res;
2911  }  }
2912    
2913  }  }

Legend:
Removed from v.757  
changed lines
  Added in v.2105

  ViewVC Help
Powered by ViewVC 1.1.26