/[escript]/trunk/escript/src/DataTagged.cpp
ViewVC logotype

Diff of /trunk/escript/src/DataTagged.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

temp_trunk_copy/escript/src/DataTagged.cpp revision 1384 by phornby, Fri Jan 11 02:29:38 2008 UTC trunk/escript/src/DataTagged.cpp revision 1946 by jfenwick, Wed Oct 29 05:48:53 2008 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2008 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15    #include "Data.h"
16  #include "DataTagged.h"  #include "DataTagged.h"
17    #include "esysUtils/esys_malloc.h"
18    
19  #include "DataConstant.h"  #include "DataConstant.h"
20  #include "DataException.h"  #include "DataException.h"
21  #ifdef USE_NETCDF  #ifdef USE_NETCDF
22  #include <netcdfcpp.h>  #include <netcdfcpp.h>
23  #endif  #endif
24    #ifdef PASO_MPI
25    #include <mpi.h>
26    #endif
27    #include "DataMaths.h"
28    
29  using namespace std;  using namespace std;
30    
31  namespace escript {  namespace escript {
32    
33  DataTagged::DataTagged()  DataTagged::DataTagged()
34    : DataAbstract(FunctionSpace())    : DataAbstract(FunctionSpace(),DataTypes::scalarShape)
35  {  {
36    // default constructor    // default constructor
37    
38    // create a scalar default value    // create a scalar default value
39    m_data.resize(1,0.,1);    m_data.resize(1,0.,1);
40    DataArrayView temp(m_data,DataArrayView::ShapeType());  /*  DataArrayView temp(m_data,DataTypes::ShapeType());
41    setPointDataView(temp);    setPointDataView(temp);*/
42  }  }
43    
44  DataTagged::DataTagged(const TagListType& tagKeys,  // DataTagged::DataTagged(const TagListType& tagKeys,
45                 const ValueListType& values,  //             const ValueListType& values,
46                 const DataArrayView& defaultValue,  //             const DataArrayView& defaultValue,
47                 const FunctionSpace& what)  //             const FunctionSpace& what)
48    : DataAbstract(what)  //   : DataAbstract(what)
49  {  // {
50    // constructor  //   // constructor
51    //
52    // initialise the array of data values  //   // initialise the array of data values
53    // the default value is always the first item in the values list  //   // the default value is always the first item in the values list
54    int len = defaultValue.noValues();  //   int len = defaultValue.noValues();
55    m_data.resize(len,0.,len);  //   m_data.resize(len,0.,len);
56    for (int i=0; i<defaultValue.noValues(); i++) {  //   for (int i=0; i<defaultValue.noValues(); i++) {
57      m_data[i]=defaultValue.getData(i);  //     m_data[i]=defaultValue.getData(i);
58    }  //   }
59    //
60    // create the data view  //   // create the data view
61    DataArrayView temp(m_data,defaultValue.getShape());  //   DataArrayView temp(m_data,defaultValue.getShape());
62    setPointDataView(temp);  //   setPointDataView(temp);
63    //
64    // add remaining tags and values  //   // add remaining tags and values
65    addTaggedValues(tagKeys,values);  //   addTaggedValues(tagKeys,values);
66  }  // }
67    
68  DataTagged::DataTagged(const FunctionSpace& what,  DataTagged::DataTagged(const FunctionSpace& what,
69                         const DataArrayView::ShapeType &shape,                         const DataTypes::ShapeType &shape,
70                         const int tags[],                         const int tags[],
71                         const ValueType& data)                         const ValueType& data)
72    : DataAbstract(what)    : DataAbstract(what,shape)
73  {  {
74    // alternative constructor    // alternative constructor
75    // not unit_tested tested yet    // not unit_tested tested yet
76      // It is not explicitly unit tested yet, but it is called from DataFactory
77    
78      if (!what.canTag())
79      {
80        throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
81      }
82    // copy the data    // copy the data
83    m_data=data;    m_data=data;
84    
85    // create the view of the data    // create the view of the data
86    DataArrayView tempView(m_data,shape);  //   DataArrayView tempView(m_data,shape);
87    setPointDataView(tempView);  //   setPointDataView(tempView);
88    
89      // we can't rely on the tag array to give us the number of tags so
90      // use the data we have been passed
91      int valsize=DataTypes::noValues(shape);
92      int ntags=data.size()/valsize;
93    
94    // create the tag lookup map    // create the tag lookup map
95    for (int sampleNo=0; sampleNo<getNumSamples(); sampleNo++) {    // we assume that the first value and first tag are the default value so we skip
96      m_offsetLookup.insert(DataMapType::value_type(sampleNo,tags[sampleNo]));    for (int i=1;i<ntags;++i)
97      {
98        m_offsetLookup.insert(DataMapType::value_type(tags[i],i*valsize));
99    }    }
100  }  }
101    
102  DataTagged::DataTagged(const FunctionSpace& what,  DataTagged::DataTagged(const FunctionSpace& what,
103                         const DataArrayView::ShapeType &shape,                         const DataTypes::ShapeType &shape,
104                         const TagListType& tags,                         const TagListType& tags,
105                         const ValueType& data)                         const ValueType& data)
106    : DataAbstract(what)    : DataAbstract(what,shape)
107  {  {
108    // alternative constructor    // alternative constructor
109    // not unit_tested tested yet  
110      if (!what.canTag())
111      {
112        throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
113      }
114    
115    // copy the data    // copy the data
116    m_data=data;    m_data=data;
117    
118    // create the view of the data    // create the view of the data
119    DataArrayView tempView(m_data,shape);  //   DataArrayView tempView(m_data,shape);
120    setPointDataView(tempView);  //   setPointDataView(tempView);
121    
122      // create the tag lookup map
123    
124    //   for (int sampleNo=0; sampleNo<getNumSamples(); sampleNo++) {
125    //     m_offsetLookup.insert(DataMapType::value_type(sampleNo,tags[sampleNo]));
126    //   }
127    
128      // The above code looks like it will create a map the wrong way around
129    
130      int valsize=DataTypes::noValues(shape);
131      int npoints=(data.size()/valsize)-1;
132      int ntags=tags.size();
133      if (ntags>npoints)
134      {     // This throw is not unit tested yet
135        throw DataException("Programming error - Too many tags for the supplied values.");
136      }
137    
138    // create the tag lookup map    // create the tag lookup map
139    for (int sampleNo=0; sampleNo<getNumSamples(); sampleNo++) {    // we assume that the first value is the default value so we skip it (hence the i+1 below)
140      m_offsetLookup.insert(DataMapType::value_type(sampleNo,tags[sampleNo]));    for (int i=0;i<ntags;++i)
141      {
142        m_offsetLookup.insert(DataMapType::value_type(tags[i],(i+1)*valsize));
143    }    }
144  }  }
145    
146    
147  DataTagged::DataTagged(const DataTagged& other)  DataTagged::DataTagged(const DataTagged& other)
148    : DataAbstract(other.getFunctionSpace()),    : DataAbstract(other.getFunctionSpace(),other.getShape()),
149    m_data(other.m_data),    m_offsetLookup(other.m_offsetLookup),
150    m_offsetLookup(other.m_offsetLookup)    m_data(other.m_data)
151  {  {
152    // copy constructor    // copy constructor
153    
154    // create the data view    // create the data view
155    DataArrayView temp(m_data,other.getPointDataView().getShape());  //   DataArrayView temp(m_data,other.getPointDataView().getShape());
156    setPointDataView(temp);  //   setPointDataView(temp);
157  }  }
158    
159  DataTagged::DataTagged(const DataConstant& other)  DataTagged::DataTagged(const DataConstant& other)
160    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace(),other.getShape())
161  {  {
162    // copy constructor    // copy constructor
163    
164      if (!other.getFunctionSpace().canTag())
165      {
166        throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
167      }
168    
169    // fill the default value with the constant value item from "other"    // fill the default value with the constant value item from "other"
170    const DataArrayView& value=other.getPointDataView();  //   const DataArrayView& value=other.getPointDataView();
171    int len = value.noValues();    int len = other.getNoValues();
172    m_data.resize(len,0.,len);    m_data.resize(len,0.,len);
173    for (int i=0; i<value.noValues(); i++) {    for (int i=0; i<len; i++) {
174      m_data[i]=value.getData(i);      m_data[i]=other.getVector()[i];
175    }    }
176    
177    // create the data view  //   // create the data view
178    DataArrayView temp(m_data,value.getShape());  //   DataArrayView temp(m_data,value.getShape());
179    setPointDataView(temp);  //   setPointDataView(temp);
180    }
181    
182    
183    // Create a new object by copying tags
184    DataTagged::DataTagged(const FunctionSpace& what,
185                 const DataTypes::ShapeType& shape,
186             const DataTypes::ValueType& defaultvalue,
187                 const DataTagged* tagsource)
188     : DataAbstract(what,shape)
189    {
190    // This constructor has not been unit tested yet
191    
192      if (defaultvalue.size()!=DataTypes::noValues(shape)) {
193        throw DataException("Programming error - defaultvalue does not match supplied shape.");
194      }
195    
196    
197      if (!what.canTag())
198      {
199        throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
200      }
201    
202      if (tagsource!=0)
203      {
204           m_data.resize(defaultvalue.size(),0.);   // since this is tagged data, we should have blocksize=1
205    
206           DataTagged::DataMapType::const_iterator i;
207           for (i=tagsource->getTagLookup().begin();i!=tagsource->getTagLookup().end();i++) {
208          addTag(i->first);
209           }
210      }
211      else
212      {
213        m_data.resize(defaultvalue.size());
214      }
215    
216      // need to set the default value ....
217      for (int i=0; i<defaultvalue.size(); i++) {
218         m_data[i]=defaultvalue[i];
219      }
220    }
221    
222    DataAbstract*
223    DataTagged::deepCopy()
224    {
225      return new DataTagged(*this);
226  }  }
227    
228  DataAbstract*  DataAbstract*
229  DataTagged::getSlice(const DataArrayView::RegionType& region) const  DataTagged::getSlice(const DataTypes::RegionType& region) const
230  {  {
231    return new DataTagged(*this, region);    return new DataTagged(*this, region);
232  }  }
233    
234  DataTagged::DataTagged(const DataTagged& other,  DataTagged::DataTagged(const DataTagged& other,
235                 const DataArrayView::RegionType& region)                 const DataTypes::RegionType& region)
236    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
237  {  {
238    // slice constructor    // slice constructor
239    
240    // get the shape of the slice to copy from other    // get the shape of the slice to copy from other
241    DataArrayView::ShapeType regionShape(DataArrayView::getResultSliceShape(region));    DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
242    DataArrayView::RegionLoopRangeType regionLoopRange=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
243    
244    // allocate enough space in this for all values    // allocate enough space in this for all values
245    // (need to add one to allow for the default value)    // (need to add one to allow for the default value)
246    int len = DataArrayView::noValues(regionShape)*(other.m_offsetLookup.size()+1);    int len = DataTypes::noValues(regionShape)*(other.m_offsetLookup.size()+1);
247    m_data.resize(len,0.0,len);    m_data.resize(len,0.0,len);
248    
249    // create the data view    // create the data view
250    DataArrayView temp(m_data,regionShape);  //   DataArrayView temp(m_data,regionShape);
251    setPointDataView(temp);  //   setPointDataView(temp);
252    
253    // copy the default value from other to this    // copy the default value from other to this
254    getDefaultValue().copySlice(other.getDefaultValue(), regionLoopRange);  //   getDefaultValue().copySlice(other.getDefaultValue(), regionLoopRange);
255      const DataTypes::ShapeType& otherShape=other.getShape();
256      const DataTypes::ValueType& otherData=other.getVector();
257      DataTypes::copySlice(getVector(),getShape(),getDefaultOffset(),otherData,otherShape,other.getDefaultOffset(), regionLoopRange);
258    
259    // loop through the tag values copying these    // loop through the tag values copying these
260    DataMapType::const_iterator pos;    DataMapType::const_iterator pos;
261    DataArrayView::ValueType::size_type tagOffset=getPointDataView().noValues();    DataTypes::ValueType::size_type tagOffset=getNoValues();
262    for (pos=other.m_offsetLookup.begin();pos!=other.m_offsetLookup.end();pos++){    for (pos=other.m_offsetLookup.begin();pos!=other.m_offsetLookup.end();pos++){
263      getPointDataView().copySlice(tagOffset,other.getPointDataView(),pos->second,regionLoopRange);  //     getPointDataView().copySlice(tagOffset,other.getPointDataView(),pos->second,regionLoopRange);
264        DataTypes::copySlice(m_data,getShape(),tagOffset,otherData, otherShape, pos->second, regionLoopRange);
265      m_offsetLookup.insert(DataMapType::value_type(pos->first,tagOffset));      m_offsetLookup.insert(DataMapType::value_type(pos->first,tagOffset));
266      tagOffset+=getPointDataView().noValues();      tagOffset+=getNoValues();
267    }    }
268  }  }
269    
270  void  void
271  DataTagged::setSlice(const DataAbstract* other,  DataTagged::setSlice(const DataAbstract* other,
272                       const DataArrayView::RegionType& region)                       const DataTypes::RegionType& region)
273  {  {
274    
275    // other must be another DataTagged object    // other must be another DataTagged object
# Line 186  DataTagged::setSlice(const DataAbstract* Line 280  DataTagged::setSlice(const DataAbstract*
280    }    }
281    
282    // determine shape of the specified region    // determine shape of the specified region
283    DataArrayView::ShapeType regionShape(DataArrayView::getResultSliceShape(region));    DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
284    
285    // modify region specification as needed to match rank of this object    // modify region specification as needed to match rank of this object
286    DataArrayView::RegionLoopRangeType regionLoopRange=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
287    
288    // ensure rank/shape of this object is compatible with specified region    // ensure rank/shape of this object is compatible with specified region
289    if (getPointDataView().getRank()!=region.size()) {    if (getRank()!=region.size()) {
290      throw DataException("Error - Invalid slice region.");      throw DataException("Error - Invalid slice region.");
291    }    }
292    if (otherTemp->getPointDataView().getRank()>0 && !other->getPointDataView().checkShape(regionShape)) {    if (otherTemp->getRank()>0 && !DataTypes::checkShape(other->getShape(),regionShape)) {
293      throw DataException (other->getPointDataView().createShapeErrorMessage(      throw DataException (DataTypes::createShapeErrorMessage(
294                           "Error - Couldn't copy slice due to shape mismatch.",regionShape));                           "Error - Couldn't copy slice due to shape mismatch.",regionShape,other->getShape()));
295    }    }
296    
297      const DataTypes::ValueType& otherData=otherTemp->getVector();
298      const DataTypes::ShapeType& otherShape=otherTemp->getShape();
299    // copy slice from other default value to this default value    // copy slice from other default value to this default value
300    getDefaultValue().copySliceFrom(otherTemp->getDefaultValue(), regionLoopRange);  //   getDefaultValue().copySliceFrom(otherTemp->getDefaultValue(), regionLoopRange);
301      DataTypes::copySliceFrom(m_data,getShape(),getDefaultOffset(),otherData,otherShape,otherTemp->getDefaultOffset(),regionLoopRange);
302    
303    // loop through tag values in other, adding any which aren't in this, using default value    // loop through tag values in other, adding any which aren't in this, using default value
304    DataMapType::const_iterator pos;    DataMapType::const_iterator pos;
305    for (pos=otherTemp->m_offsetLookup.begin();pos!=otherTemp->m_offsetLookup.end();pos++) {    for (pos=otherTemp->m_offsetLookup.begin();pos!=otherTemp->m_offsetLookup.end();pos++) {
306      if (!isCurrentTag(pos->first)) {      if (!isCurrentTag(pos->first)) {
307        addTaggedValue(pos->first,getDefaultValue());        addTag(pos->first);
308      }      }
309    }    }
310    
311    // loop through the tag values copying slices from other to this    // loop through the tag values copying slices from other to this
312    for (pos=m_offsetLookup.begin();pos!=m_offsetLookup.end();pos++) {    for (pos=m_offsetLookup.begin();pos!=m_offsetLookup.end();pos++) {
313      getDataPointByTag(pos->first).copySliceFrom(otherTemp->getDataPointByTag(pos->first), regionLoopRange);  //     getDataPointByTag(pos->first).copySliceFrom(otherTemp->getDataPointByTag(pos->first), regionLoopRange);
314        DataTypes::copySliceFrom(m_data,getShape(),getOffsetForTag(pos->first),otherData, otherShape, otherTemp->getOffsetForTag(pos->first), regionLoopRange);
315    
316    }    }
317    
318  }  }
# Line 248  DataTagged::getTagNumber(int dpno) Line 347  DataTagged::getTagNumber(int dpno)
347    return(tagNo);    return(tagNo);
348  }  }
349    
350    // void
351    // DataTagged::setTaggedValues(const TagListType& tagKeys,
352    //                             const ValueListType& values)
353    // {
354    //   addTaggedValues(tagKeys,values);
355    // }
356    
357  void  void
358  DataTagged::setTaggedValues(const TagListType& tagKeys,  DataTagged::setTaggedValue(int tagKey,
359                              const ValueListType& values)                 const DataTypes::ShapeType& pointshape,
360  {                             const ValueType& value,
361    addTaggedValues(tagKeys,values);                 int dataOffset)
362    {
363      if (!DataTypes::checkShape(getShape(), pointshape)) {
364          throw DataException(DataTypes::createShapeErrorMessage(
365                              "Error - Cannot setTaggedValue due to shape mismatch.", pointshape,getShape()));
366      }
367      DataMapType::iterator pos(m_offsetLookup.find(tagKey));
368      if (pos==m_offsetLookup.end()) {
369        // tag couldn't be found so use addTaggedValue
370        addTaggedValue(tagKey,pointshape, value, dataOffset);
371      } else {
372        // copy the values into the data array at the offset determined by m_offsetLookup
373        int offset=pos->second;
374        for (int i=0; i<getNoValues(); i++) {
375          m_data[offset+i]=value[i+dataOffset];
376        }
377      }
378  }  }
379    
380    
381    /*
382  void  void
383  DataTagged::setTaggedValue(int tagKey,  DataTagged::setTaggedValue(int tagKey,
384                             const DataArrayView& value)                             const DataArrayView& value)
385  {  {
386    if (!getPointDataView().checkShape(value.getShape())) {    if (!getPointDataView().checkShape(value.getShape())) {
387        throw DataException(getPointDataView().createShapeErrorMessage(        throw DataException(DataTypes::createShapeErrorMessage(
388                            "Error - Cannot setTaggedValue due to shape mismatch.", value.getShape()));                            "Error - Cannot setTaggedValue due to shape mismatch.", value.getShape(),getShape()));
389    }    }
390    DataMapType::iterator pos(m_offsetLookup.find(tagKey));    DataMapType::iterator pos(m_offsetLookup.find(tagKey));
391    if (pos==m_offsetLookup.end()) {    if (pos==m_offsetLookup.end()) {
# Line 274  DataTagged::setTaggedValue(int tagKey, Line 398  DataTagged::setTaggedValue(int tagKey,
398        m_data[offset+i]=value.getData(i);        m_data[offset+i]=value.getData(i);
399      }      }
400    }    }
401    }*/
402    
403    // void
404    // DataTagged::addTaggedValues(const TagListType& tagKeys,
405    //                             const ValueListType& values)
406    // {
407    //   if (values.size()==0) {
408    //     // copy the current default value for each of the tags
409    //     TagListType::const_iterator iT;
410    //     for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
411    //       // the point data view for DataTagged points at the default value
412    //       addTaggedValue(*iT,getPointDataView());
413    //     }
414    //   } else if (values.size()==1 && tagKeys.size()>1) {
415    //     // assume the one given value will be used for all tag values
416    //     TagListType::const_iterator iT;
417    //     for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
418    //       addTaggedValue(*iT,values[0]);
419    //     }
420    //   } else {
421    //     if (tagKeys.size()!=values.size()) {
422    //       stringstream temp;
423    //       temp << "Error - (addTaggedValue) Number of tags: " << tagKeys.size()
424    //     << " doesn't match number of values: " << values.size();
425    //       throw DataException(temp.str());
426    //     } else {
427    //       unsigned int i;
428    //       for (i=0;i<tagKeys.size();i++) {
429    //         addTaggedValue(tagKeys[i],values[i]);
430    //       }
431    //     }
432    //   }
433    // }
434    
435    
436    void
437    DataTagged::addTaggedValues(const TagListType& tagKeys,
438                                const ValueBatchType& values,
439                                const ShapeType& vShape)
440    {
441      DataTypes::ValueType t(values.size(),0);
442      for (size_t i=0;i<values.size();++i)
443      {
444        t[i]=values[i];
445      }
446      addTaggedValues(tagKeys,t,vShape);
447  }  }
448    
449    
450    // Note: The check to see if vShape==our shape is done in the addTaggedValue method
451  void  void
452  DataTagged::addTaggedValues(const TagListType& tagKeys,  DataTagged::addTaggedValues(const TagListType& tagKeys,
453                              const ValueListType& values)                              const ValueType& values,
454                                const ShapeType& vShape)
455  {  {
456      int n=getNoValues();
457      int numVals=values.size()/n;
458    if (values.size()==0) {    if (values.size()==0) {
459      // copy the current default value for each of the tags      // copy the current default value for each of the tags
460      TagListType::const_iterator iT;      TagListType::const_iterator iT;
461      for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {      for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
462        // the point data view for DataTagged points at the default value        // the point data view for DataTagged points at the default value
463        addTaggedValue(*iT,getPointDataView());        addTag(*iT);
464      }      }
465    } else if (values.size()==1 && tagKeys.size()>1) {    } else if (numVals==1 && tagKeys.size()>1) {
466      // assume the one given value will be used for all tag values      // assume the one given value will be used for all tag values
467      TagListType::const_iterator iT;      TagListType::const_iterator iT;
468      for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {      for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
469        addTaggedValue(*iT,values[0]);        addTaggedValue(*iT, vShape, values,0);
470      }      }
471    } else {    } else {
472      if (tagKeys.size()!=values.size()) {      if (tagKeys.size()!=numVals) {
473        stringstream temp;        stringstream temp;
474        temp << "Error - (addTaggedValue) Number of tags: " << tagKeys.size()        temp << "Error - (addTaggedValue) Number of tags: " << tagKeys.size()
475         << " doesn't match number of values: " << values.size();         << " doesn't match number of values: " << values.size();
476        throw DataException(temp.str());        throw DataException(temp.str());
477      } else {      } else {
478        for (int i=0;i<tagKeys.size();i++) {        unsigned int i;
479          addTaggedValue(tagKeys[i],values[i]);        int offset=0;
480          for (i=0;i<tagKeys.size();i++ ,offset+=n) {
481            addTaggedValue(tagKeys[i],vShape,values,offset);
482        }        }
483      }      }
484    }    }
485  }  }
486    
487    
488    
489    
490  void  void
491  DataTagged::addTaggedValue(int tagKey,  DataTagged::addTaggedValue(int tagKey,
492                             const DataArrayView& value)                 const DataTypes::ShapeType& pointshape,
493  {                             const ValueType& value,
494    if (!getPointDataView().checkShape(value.getShape())) {                 int dataOffset)
495      throw DataException(getPointDataView().createShapeErrorMessage(  {
496                          "Error - Cannot addTaggedValue due to shape mismatch.", value.getShape()));    if (!DataTypes::checkShape(getShape(), pointshape)) {
497        throw DataException(DataTypes::createShapeErrorMessage(
498                            "Error - Cannot addTaggedValue due to shape mismatch.", pointshape,getShape()));
499    }    }
500    DataMapType::iterator pos(m_offsetLookup.find(tagKey));    DataMapType::iterator pos(m_offsetLookup.find(tagKey));
501    if (pos!=m_offsetLookup.end()) {    if (pos!=m_offsetLookup.end()) {
502      // tag already exists so use setTaggedValue      // tag already exists so use setTaggedValue
503      setTaggedValue(tagKey,value);      setTaggedValue(tagKey,pointshape, value, dataOffset);
504    } else {    } else {
505      // save the key and the location of its data in the lookup tab      // save the key and the location of its data in the lookup tab
506      m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));      m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));
# Line 327  DataTagged::addTaggedValue(int tagKey, Line 509  DataTagged::addTaggedValue(int tagKey,
509      // all the old values plus the value to be added back into m_data      // all the old values plus the value to be added back into m_data
510      ValueType m_data_temp(m_data);      ValueType m_data_temp(m_data);
511      int oldSize=m_data.size();      int oldSize=m_data.size();
512      int newSize=m_data.size()+value.noValues();      int newSize=m_data.size()+getNoValues();
513      m_data.resize(newSize,0.,newSize);      m_data.resize(newSize,0.,newSize);
514      for (int i=0;i<oldSize;i++) {      for (int i=0;i<oldSize;i++) {
515        m_data[i]=m_data_temp[i];        m_data[i]=m_data_temp[i];
516      }      }
517      for (int i=0;i<value.noValues();i++) {      for (int i=0;i<getNoValues();i++) {
518        m_data[oldSize+i]=value.getData(i);        m_data[oldSize+i]=value[i+dataOffset];
519      }      }
520    }    }
521  }  }
522    
523    
524    
525    
526    // void
527    // DataTagged::addTaggedValue(int tagKey,
528    //                            const DataArrayView& value)
529    // {
530    //   if (!getPointDataView().checkShape(value.getShape())) {
531    //     throw DataException(DataTypes::createShapeErrorMessage(
532    //                         "Error - Cannot addTaggedValue due to shape mismatch.", value.getShape(),getShape()));
533    //   }
534    //   DataMapType::iterator pos(m_offsetLookup.find(tagKey));
535    //   if (pos!=m_offsetLookup.end()) {
536    //     // tag already exists so use setTaggedValue
537    //     setTaggedValue(tagKey,value);
538    //   } else {
539    //     // save the key and the location of its data in the lookup tab
540    //     m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));
541    //     // add the data given in "value" at the end of m_data
542    //     // need to make a temp copy of m_data, resize m_data, then copy
543    //     // all the old values plus the value to be added back into m_data
544    //     ValueType m_data_temp(m_data);
545    //     int oldSize=m_data.size();
546    //     int newSize=m_data.size()+value.noValues();
547    //     m_data.resize(newSize,0.,newSize);
548    //     for (int i=0;i<oldSize;i++) {
549    //       m_data[i]=m_data_temp[i];
550    //     }
551    //     for (int i=0;i<value.noValues();i++) {
552    //       m_data[oldSize+i]=value.getData(i);
553    //     }
554    //   }
555    // }
556    
557    
558    void
559    DataTagged::addTag(int tagKey)
560    {
561      DataMapType::iterator pos(m_offsetLookup.find(tagKey));
562      if (pos!=m_offsetLookup.end()) {
563        // tag already exists so use setTaggedValue
564    //    setTaggedValue(tagKey,value);
565      } else {
566        // save the key and the location of its data in the lookup tab
567        m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));
568        // add the data given in "value" at the end of m_data
569        // need to make a temp copy of m_data, resize m_data, then copy
570        // all the old values plus the value to be added back into m_data
571        ValueType m_data_temp(m_data);
572        int oldSize=m_data.size();
573        int newSize=m_data.size()+getNoValues();
574        m_data.resize(newSize,0.,newSize);
575        for (int i=0;i<oldSize;i++) {
576          m_data[i]=m_data_temp[i];
577        }
578        for (int i=0;i<getNoValues();i++) {
579          m_data[oldSize+i]=m_data[m_defaultValueOffset+i];
580        }
581      }
582    }
583    
584    
585  double*  double*
586  DataTagged::getSampleDataByTag(int tag)  DataTagged::getSampleDataByTag(int tag)
587  {  {
# Line 354  DataTagged::getSampleDataByTag(int tag) Line 598  DataTagged::getSampleDataByTag(int tag)
598  string  string
599  DataTagged::toString() const  DataTagged::toString() const
600  {  {
601      using namespace escript::DataTypes;
602      string empty="";
603    stringstream temp;    stringstream temp;
604    DataMapType::const_iterator i;    DataMapType::const_iterator i;
605    temp << "Tag(Default)" << endl;    temp << "Tag(Default)" << endl;
606    temp << getDefaultValue().toString() << endl;    temp << pointToString(m_data,getShape(),getDefaultOffset(),empty) << endl;
607    // create a temporary view as the offset will be changed    // create a temporary view as the offset will be changed
608    DataArrayView tempView(getPointDataView().getData(), getPointDataView().getShape());  //   DataArrayView tempView(getPointDataView().getData(), getPointDataView().getShape());
609    for (i=m_offsetLookup.begin();i!=m_offsetLookup.end();++i) {    for (i=m_offsetLookup.begin();i!=m_offsetLookup.end();++i) {
610      temp << "Tag(" << i->first << ")" << endl;      temp << "Tag(" << i->first << ")" << endl;
611      tempView.setOffset(i->second);      temp << pointToString(m_data,getShape(),i->second,empty) << endl;
612      temp << tempView.toString() << endl;  //     tempView.setOffset(i->second);
613    //     temp << tempView.toString() << endl;
614    }    }
615    return temp.str();    return temp.str();
616  }  }
617    
618  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
619  DataTagged::getPointOffset(int sampleNo,  DataTagged::getPointOffset(int sampleNo,
620                             int dataPointNo) const                             int dataPointNo) const
621  {  {
622    int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);    int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
623    DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));    DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));
624    DataArrayView::ValueType::size_type offset=m_defaultValueOffset;    DataTypes::ValueType::size_type offset=m_defaultValueOffset;
625    if (pos!=m_offsetLookup.end()) {    if (pos!=m_offsetLookup.end()) {
626      offset=pos->second;      offset=pos->second;
627    }    }
628    return offset;    return offset;
629  }  }
630    
631  DataArrayView  // DataArrayView
632  DataTagged::getDataPointByTag(int tag) const  // DataTagged::getDataPointByTag(int tag) const
633    // {
634    //   DataMapType::const_iterator pos(m_offsetLookup.find(tag));
635    //   DataTypes::ValueType::size_type offset=m_defaultValueOffset;
636    //   if (pos!=m_offsetLookup.end()) {
637    //     offset=pos->second;
638    //   }
639    //   DataArrayView temp(getPointDataView());
640    //   temp.setOffset(offset);
641    //   return temp;
642    // }
643    //
644    
645    
646    DataTypes::ValueType::size_type
647    DataTagged::getOffsetForTag(int tag) const
648  {  {
649    DataMapType::const_iterator pos(m_offsetLookup.find(tag));    DataMapType::const_iterator pos(m_offsetLookup.find(tag));
650    DataArrayView::ValueType::size_type offset=m_defaultValueOffset;    DataTypes::ValueType::size_type offset=m_defaultValueOffset;
651    if (pos!=m_offsetLookup.end()) {    if (pos!=m_offsetLookup.end()) {
652      offset=pos->second;      offset=pos->second;
653    }    }
654    DataArrayView temp(getPointDataView());    return offset;
   temp.setOffset(offset);  
   return temp;  
655  }  }
656    
657  DataArrayView  DataTypes::ValueType::const_reference
658  DataTagged::getDataPoint(int sampleNo,  DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i) const
                          int dataPointNo)  
659  {  {
660    EsysAssert(validSampleNo(sampleNo),"(getDataPoint) Invalid sampleNo: " << sampleNo);    DataMapType::const_iterator pos(m_offsetLookup.find(tag));
661    int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);    DataTypes::ValueType::size_type offset=m_defaultValueOffset;
662    return getDataPointByTag(tagKey);    if (pos!=m_offsetLookup.end()) {
663        offset=pos->second;
664      }
665      return m_data[offset+i];
666    /*  DataArrayView temp(getPointDataView());
667      temp.setOffset(offset);
668      return temp.getData()[offset+i];*/
669  }  }
670    
 int  
 DataTagged::archiveData(ofstream& archiveFile,  
                         const DataArrayView::ValueType::size_type noValues) const  
 {  
   return(m_data.archiveData(archiveFile, noValues));  
 }  
671    
672  int  DataTypes::ValueType::reference
673  DataTagged::extractData(ifstream& archiveFile,  DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i)
                         const DataArrayView::ValueType::size_type noValues)  
674  {  {
675    return(m_data.extractData(archiveFile, noValues));    DataMapType::const_iterator pos(m_offsetLookup.find(tag));
676      DataTypes::ValueType::size_type offset=m_defaultValueOffset;
677      if (pos!=m_offsetLookup.end()) {
678        offset=pos->second;
679      }
680      return m_data[offset+i];
681    /*  DataArrayView temp(getPointDataView());
682      temp.setOffset(offset);
683      return temp.getData()[offset+i];*/
684  }  }
685    
686    
687    
688    
689    
690    
691    // DataArrayView
692    // DataTagged::getDataPoint(int sampleNo,
693    //                          int dataPointNo)
694    // {
695    //   EsysAssert(validSampleNo(sampleNo),"(getDataPoint) Invalid sampleNo: " << sampleNo);
696    //   int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
697    //   return getDataPointByTag(tagKey);
698    // }
699    
700    
701  void  void
702  DataTagged::symmetric(DataAbstract* ev)  DataTagged::symmetric(DataAbstract* ev)
703  {  {
# Line 426  DataTagged::symmetric(DataAbstract* ev) Line 708  DataTagged::symmetric(DataAbstract* ev)
708    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
709    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
710    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
711      ValueType& evVec=temp_ev->getVector();
712      const ShapeType& evShape=temp_ev->getShape();
713    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
714        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
715        DataArrayView thisView=getDataPointByTag(i->first);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
716        DataArrayView evView=temp_ev->getDataPointByTag(i->first);  //       DataArrayView thisView=getDataPointByTag(i->first);
717        DataArrayView::symmetric(thisView,0,evView,0);  //       DataArrayView evView=temp_ev->getDataPointByTag(i->first);
718          DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
719    
720    //       DataArrayView::symmetric(thisView,0,evView,0);
721          DataMaths::symmetric(m_data,getShape(),offset,evVec, evShape, evoffset);
722    }    }
723    DataArrayView::symmetric(getDefaultValue(),0,temp_ev->getDefaultValue(),0);  //   symmetric(m_data,getShape(),getDefaultOffset(),
724      DataMaths::symmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
725  }  }
726    
727    
728  void  void
729  DataTagged::nonsymmetric(DataAbstract* ev)  DataTagged::nonsymmetric(DataAbstract* ev)
730  {  {
# Line 444  DataTagged::nonsymmetric(DataAbstract* e Line 735  DataTagged::nonsymmetric(DataAbstract* e
735    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
736    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
737    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
738      ValueType& evVec=temp_ev->getVector();
739      const ShapeType& evShape=temp_ev->getShape();
740    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
741        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
742        DataArrayView thisView=getDataPointByTag(i->first);  /*      DataArrayView thisView=getDataPointByTag(i->first);
743        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataArrayView evView=temp_ev->getDataPointByTag(i->first);*/
744        DataArrayView::nonsymmetric(thisView,0,evView,0);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
745          DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
746          DataMaths::nonsymmetric(m_data,getShape(),offset,evVec, evShape, evoffset);
747    }    }
748    DataArrayView::nonsymmetric(getDefaultValue(),0,temp_ev->getDefaultValue(),0);    DataMaths::nonsymmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
749  }  }
750    
751    
752  void  void
753  DataTagged::trace(DataAbstract* ev, int axis_offset)  DataTagged::trace(DataAbstract* ev, int axis_offset)
754  {  {
# Line 462  DataTagged::trace(DataAbstract* ev, int Line 759  DataTagged::trace(DataAbstract* ev, int
759    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
760    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
761    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
762      ValueType& evVec=temp_ev->getVector();
763      const ShapeType& evShape=temp_ev->getShape();
764    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
765        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
766        DataArrayView thisView=getDataPointByTag(i->first);  //       DataArrayView thisView=getDataPointByTag(i->first);
767        DataArrayView evView=temp_ev->getDataPointByTag(i->first);  //       DataArrayView evView=temp_ev->getDataPointByTag(i->first);
768        DataArrayView::trace(thisView,0,evView,0, axis_offset);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
769          DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
770          DataMaths::trace(m_data,getShape(),offset,evVec, evShape, evoffset, axis_offset);
771    }    }
772    DataArrayView::trace(getDefaultValue(),0,temp_ev->getDefaultValue(),0,axis_offset);    DataMaths::trace(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
773  }  }
774    
775  void  void
# Line 481  DataTagged::transpose(DataAbstract* ev, Line 782  DataTagged::transpose(DataAbstract* ev,
782    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
783    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
784    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
785      ValueType& evVec=temp_ev->getVector();
786      const ShapeType& evShape=temp_ev->getShape();
787    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
788        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
789        DataArrayView thisView=getDataPointByTag(i->first);  //       DataArrayView thisView=getDataPointByTag(i->first);
790        DataArrayView evView=temp_ev->getDataPointByTag(i->first);  //       DataArrayView evView=temp_ev->getDataPointByTag(i->first);
791        DataArrayView::transpose(thisView,0,evView,0, axis_offset);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
792          DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
793          DataMaths::transpose(m_data,getShape(),offset,evVec, evShape, evoffset, axis_offset);
794    }    }
795    DataArrayView::transpose(getDefaultValue(),0,temp_ev->getDefaultValue(),0,axis_offset);    DataMaths::transpose(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
796  }  }
797    
798  void  void
# Line 500  DataTagged::swapaxes(DataAbstract* ev, i Line 805  DataTagged::swapaxes(DataAbstract* ev, i
805    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
806    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
807    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
808      ValueType& evVec=temp_ev->getVector();
809      const ShapeType& evShape=temp_ev->getShape();
810    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
811        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
812        DataArrayView thisView=getDataPointByTag(i->first);  /*      DataArrayView thisView=getDataPointByTag(i->first);
813        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataArrayView evView=temp_ev->getDataPointByTag(i->first);*/
814        DataArrayView::swapaxes(thisView,0,evView,0,axis0,axis1);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
815          DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
816          DataMaths::swapaxes(m_data,getShape(),offset,evVec, evShape, evoffset,axis0,axis1);
817    }    }
818    DataArrayView::swapaxes(getDefaultValue(),0,temp_ev->getDefaultValue(),0,axis0,axis1);    DataMaths::swapaxes(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis0,axis1);
819  }  }
820    
821  void  void
# Line 519  DataTagged::eigenvalues(DataAbstract* ev Line 828  DataTagged::eigenvalues(DataAbstract* ev
828    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
829    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
830    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
831      ValueType& evVec=temp_ev->getVector();
832      const ShapeType& evShape=temp_ev->getShape();
833    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
834        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
835        DataArrayView thisView=getDataPointByTag(i->first);  //       DataArrayView thisView=getDataPointByTag(i->first);
836        DataArrayView evView=temp_ev->getDataPointByTag(i->first);  //       DataArrayView evView=temp_ev->getDataPointByTag(i->first);
837        DataArrayView::eigenvalues(thisView,0,evView,0);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
838          DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
839          DataMaths::eigenvalues(m_data,getShape(),offset,evVec, evShape, evoffset);
840    }    }
841    DataArrayView::eigenvalues(getDefaultValue(),0,temp_ev->getDefaultValue(),0);    DataMaths::eigenvalues(m_data,getShape(),getDefaultOffset(),evVec, evShape, temp_ev->getDefaultOffset());
842  }  }
843  void  void
844  DataTagged::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)  DataTagged::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
# Line 541  DataTagged::eigenvalues_and_eigenvectors Line 854  DataTagged::eigenvalues_and_eigenvectors
854    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
855    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
856    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
857      ValueType& evVec=temp_ev->getVector();
858      const ShapeType& evShape=temp_ev->getShape();
859      ValueType& VVec=temp_V->getVector();
860      const ShapeType& VShape=temp_V->getShape();
861    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
862        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
863        temp_V->addTaggedValue(i->first,temp_V->getDefaultValue());        temp_V->addTag(i->first);
864        DataArrayView thisView=getDataPointByTag(i->first);  /*      DataArrayView thisView=getDataPointByTag(i->first);
865        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataArrayView evView=temp_ev->getDataPointByTag(i->first);
866        DataArrayView VView=temp_V->getDataPointByTag(i->first);        DataArrayView VView=temp_V->getDataPointByTag(i->first);*/
867        DataArrayView::eigenvalues_and_eigenvectors(thisView,0,evView,0,VView,0,tol);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
868          DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
869          DataTypes::ValueType::size_type Voffset=temp_V->getOffsetForTag(i->first);
870    /*      DataArrayView::eigenvalues_and_eigenvectors(thisView,0,evView,0,VView,0,tol);*/
871          DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),offset,evVec, evShape, evoffset,VVec,VShape,Voffset,tol);
872    
873    }    }
874    DataArrayView::eigenvalues_and_eigenvectors(getDefaultValue(),0,    DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),getDefaultOffset(),evVec, evShape,
875                                                temp_ev->getDefaultValue(),0,                        temp_ev->getDefaultOffset(),VVec,VShape,
876                                                temp_V->getDefaultValue(),0,                        temp_V->getDefaultOffset(), tol);
                                               tol);  
877    
878    
879  }  }
880    
881  void  void
882  DataTagged::setToZero(){  DataTagged::setToZero(){
883      DataArrayView::ValueType::size_type n=m_data.size();      DataTypes::ValueType::size_type n=m_data.size();
884      for (int i=0; i<n ;++i) m_data[i]=0.;      for (int i=0; i<n ;++i) m_data[i]=0.;
885  }  }
886    
887  void  void
888  DataTagged::dump(const std::string fileName) const  DataTagged::dump(const std::string fileName) const
889  {  {
    #ifdef PASO_MPI  
    throw DataException("Error - DataTagged:: dump is not implemented for MPI yet.");  
    #endif  
890     #ifdef USE_NETCDF     #ifdef USE_NETCDF
891     const int ldims=DataArrayView::maxRank+1;     const int ldims=DataTypes::maxRank+1;
892     const NcDim* ncdims[ldims];     const NcDim* ncdims[ldims];
893     NcVar *var, *tags_var;     NcVar *var, *tags_var;
894     int rank = getPointDataView().getRank();     int rank = getRank();
895     int type=  getFunctionSpace().getTypeCode();     int type=  getFunctionSpace().getTypeCode();
896     int ndims =0;     int ndims =0;
897     long dims[ldims];     long dims[ldims];
898     const double* d_ptr=&(m_data[0]);     const double* d_ptr=&(m_data[0]);
899     DataArrayView::ShapeType shape = getPointDataView().getShape();     DataTypes::ShapeType shape = getShape();
900       int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
901       int mpi_num=getFunctionSpace().getDomain()->getMPISize();
902    #ifdef PASO_MPI
903       MPI_Status status;
904    #endif
905    
906    #ifdef PASO_MPI
907       /* Serialize NetCDF I/O */
908       if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81803, MPI_COMM_WORLD, &status);
909    #endif
910    
911     // netCDF error handler     // netCDF error handler
912     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
913     // Create the file.     // Create the file.
914     NcFile dataFile(fileName.c_str(), NcFile::Replace);     char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
915       NcFile dataFile(newFileName, NcFile::Replace);
916     // check if writing was successful     // check if writing was successful
917     if (!dataFile.is_valid())     if (!dataFile.is_valid())
918          throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");          throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");
# Line 597  DataTagged::dump(const std::string fileN Line 926  DataTagged::dump(const std::string fileN
926     if ( rank >0 ) {     if ( rank >0 ) {
927         dims[0]=shape[0];         dims[0]=shape[0];
928         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
929              throw DataException("Error - DataTagged:: appending ncdimsion 0 to netCDF file failed.");              throw DataException("Error - DataTagged:: appending ncdimension 0 to netCDF file failed.");
930     }     }
931     if ( rank >1 ) {     if ( rank >1 ) {
932         dims[1]=shape[1];         dims[1]=shape[1];
933         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
934              throw DataException("Error - DataTagged:: appending ncdimsion 1 to netCDF file failed.");              throw DataException("Error - DataTagged:: appending ncdimension 1 to netCDF file failed.");
935     }     }
936     if ( rank >2 ) {     if ( rank >2 ) {
937         dims[2]=shape[2];         dims[2]=shape[2];
938         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
939              throw DataException("Error - DataTagged:: appending ncdimsion 2 to netCDF file failed.");              throw DataException("Error - DataTagged:: appending ncdimension 2 to netCDF file failed.");
940     }     }
941     if ( rank >3 ) {     if ( rank >3 ) {
942         dims[3]=shape[3];         dims[3]=shape[3];
943         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
944              throw DataException("Error - DataTagged:: appending ncdimsion 3 to netCDF file failed.");              throw DataException("Error - DataTagged:: appending ncdimension 3 to netCDF file failed.");
945     }     }
946     const DataTagged::DataMapType& thisLookup=getTagLookup();     const DataTagged::DataMapType& thisLookup=getTagLookup();
947     DataTagged::DataMapType::const_iterator i;     DataTagged::DataMapType::const_iterator i;
948     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
949     int ntags=1;     int ntags=1;
950     for (i=thisLookup.begin();i!=thisLookupEnd;i++) ntags++;     for (i=thisLookup.begin();i!=thisLookupEnd;i++) ntags++;
951     int* tags =(int*) malloc(ntags*sizeof(int));     int* tags =(int*) esysUtils::malloc(ntags*sizeof(int));
952     int c=1;     int c=1;
953     tags[0]=-1;     tags[0]=-1;
954     for (i=thisLookup.begin();i!=thisLookupEnd;i++) tags[c++]=i->first;     for (i=thisLookup.begin();i!=thisLookupEnd;i++) tags[c++]=i->first;
955     dims[rank]=ntags;     dims[rank]=ntags;
956     if (! (ncdims[rank] = dataFile.add_dim("num_tags", dims[rank])) )     if (! (ncdims[rank] = dataFile.add_dim("num_tags", dims[rank])) )
957     {     {
958         free(tags);         esysUtils::free(tags);
959             throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");             throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");
960     }     }
961     if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )     if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )
962     {     {
963      free(tags);      esysUtils::free(tags);
964          throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");          throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");
965     }     }
966     if (! (tags_var->put(tags,dims[rank])) )     if (! (tags_var->put(tags,dims[rank])) )
967     {     {
968      free(tags);      esysUtils::free(tags);
969          throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");          throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");
970     }     }
971     if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )     if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
972     {     {
973      free(tags);      esysUtils::free(tags);
974          throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");          throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");
975     }     }
976     if (! (var->put(d_ptr,dims)) )     if (! (var->put(d_ptr,dims)) )
977     {     {
978      free(tags);      esysUtils::free(tags);
979          throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");          throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");
980     }     }
981    #ifdef PASO_MPI
982       if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81803, MPI_COMM_WORLD);
983    #endif
984     #else     #else
985     throw DataException("Error - DataTagged:: dump is not configured with netCDF. Please contact your installation manager.");     throw DataException("Error - DataTagged:: dump is not configured with netCDF. Please contact your installation manager.");
986     #endif     #endif
987  }  }
988    
989    DataTypes::ValueType&
990    DataTagged::getVector()
991    {
992        return m_data;
993    }
994    
995    const DataTypes::ValueType&
996    DataTagged::getVector() const
997    {
998        return m_data;
999    }
1000    
1001  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1384  
changed lines
  Added in v.1946

  ViewVC Help
Powered by ViewVC 1.1.26