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

revision 1023 by gross, Tue Mar 13 05:43:36 2007 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1 
 // $Id$  
1    
2  /*  /*******************************************************
3   ************************************************************  *
4   *          Copyright 2006 by ACcESS MNRF                   *  * Copyright (c) 2003-2009 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    #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    
30    #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {throw DataException("Attempt to modify shared object");}
31    
32    // #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {std::ostringstream ss; ss << " Attempt to modify shared object. line " << __LINE__ << " of " << __FILE__; throw DataException(ss.str());}
33    
34  using namespace std;  using namespace std;
35    
36  namespace escript {  namespace escript {
37    
38  DataTagged::DataTagged()  DataTagged::DataTagged()
39    : DataAbstract(FunctionSpace())    : parent(FunctionSpace(),DataTypes::scalarShape)
40  {  {
41    // default constructor    // default constructor
42    
43    // create a scalar default value    // create a scalar default value
44    m_data.resize(1,0.,1);    m_data.resize(1,0.,1);
45    DataArrayView temp(m_data,DataArrayView::ShapeType());  /*  DataArrayView temp(m_data,DataTypes::ShapeType());
46    setPointDataView(temp);    setPointDataView(temp);*/
47  }  }
48    
49  DataTagged::DataTagged(const TagListType& tagKeys,  // DataTagged::DataTagged(const TagListType& tagKeys,
50                 const ValueListType& values,  //             const ValueListType& values,
51                 const DataArrayView& defaultValue,  //             const DataArrayView& defaultValue,
52                 const FunctionSpace& what)  //             const FunctionSpace& what)
53    : DataAbstract(what)  //   : DataAbstract(what)
54  {  // {
55    // constructor  //   // constructor
56    //
57    // initialise the array of data values  //   // initialise the array of data values
58    // the default value is always the first item in the values list  //   // the default value is always the first item in the values list
59    int len = defaultValue.noValues();  //   int len = defaultValue.noValues();
60    m_data.resize(len,0.,len);  //   m_data.resize(len,0.,len);
61    for (int i=0; i<defaultValue.noValues(); i++) {  //   for (int i=0; i<defaultValue.noValues(); i++) {
62      m_data[i]=defaultValue.getData(i);  //     m_data[i]=defaultValue.getData(i);
63    }  //   }
64    //
65    // create the data view  //   // create the data view
66    DataArrayView temp(m_data,defaultValue.getShape());  //   DataArrayView temp(m_data,defaultValue.getShape());
67    setPointDataView(temp);  //   setPointDataView(temp);
68    //
69    // add remaining tags and values  //   // add remaining tags and values
70    addTaggedValues(tagKeys,values);  //   addTaggedValues(tagKeys,values);
71  }  // }
72    
73  DataTagged::DataTagged(const FunctionSpace& what,  DataTagged::DataTagged(const FunctionSpace& what,
74                         const DataArrayView::ShapeType &shape,                         const DataTypes::ShapeType &shape,
75                         const int tags[],                         const int tags[],
76                         const ValueType& data)                         const ValueType& data)
77    : DataAbstract(what)    : parent(what,shape)
78  {  {
79    // alternative constructor    // alternative constructor
80    // not unit_tested tested yet    // not unit_tested tested yet
81      // It is not explicitly unit tested yet, but it is called from DataFactory
82    
83      if (!what.canTag())
84      {
85        throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
86      }
87    // copy the data    // copy the data
88    m_data=data;    m_data=data;
89    
90    // create the view of the data    // create the view of the data
91    DataArrayView tempView(m_data,shape);  //   DataArrayView tempView(m_data,shape);
92    setPointDataView(tempView);  //   setPointDataView(tempView);
93    
94      // we can't rely on the tag array to give us the number of tags so
95      // use the data we have been passed
96      int valsize=DataTypes::noValues(shape);
97      int ntags=data.size()/valsize;
98    
99    // create the tag lookup map    // create the tag lookup map
100    for (int sampleNo=0; sampleNo<getNumSamples(); sampleNo++) {    // we assume that the first value and first tag are the default value so we skip
101      m_offsetLookup.insert(DataMapType::value_type(sampleNo,tags[sampleNo]));    for (int i=1;i<ntags;++i)
102      {
103        m_offsetLookup.insert(DataMapType::value_type(tags[i],i*valsize));
104    }    }
105  }  }
106    
107  DataTagged::DataTagged(const FunctionSpace& what,  DataTagged::DataTagged(const FunctionSpace& what,
108                         const DataArrayView::ShapeType &shape,                         const DataTypes::ShapeType &shape,
109                         const TagListType& tags,                         const TagListType& tags,
110                         const ValueType& data)                         const ValueType& data)
111    : DataAbstract(what)    : parent(what,shape)
112  {  {
113    // alternative constructor    // alternative constructor
114    // not unit_tested tested yet  
115      if (!what.canTag())
116      {
117        throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
118      }
119    
120    // copy the data    // copy the data
121    m_data=data;    m_data=data;
122    
123    // create the view of the data    // create the view of the data
124    DataArrayView tempView(m_data,shape);  //   DataArrayView tempView(m_data,shape);
125    setPointDataView(tempView);  //   setPointDataView(tempView);
126    
127      // create the tag lookup map
128    
129    //   for (int sampleNo=0; sampleNo<getNumSamples(); sampleNo++) {
130    //     m_offsetLookup.insert(DataMapType::value_type(sampleNo,tags[sampleNo]));
131    //   }
132    
133      // The above code looks like it will create a map the wrong way around
134    
135      int valsize=DataTypes::noValues(shape);
136      int npoints=(data.size()/valsize)-1;
137      int ntags=tags.size();
138      if (ntags>npoints)
139      {     // This throw is not unit tested yet
140        throw DataException("Programming error - Too many tags for the supplied values.");
141      }
142    
143    // create the tag lookup map    // create the tag lookup map
144    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)
145      m_offsetLookup.insert(DataMapType::value_type(sampleNo,tags[sampleNo]));    for (int i=0;i<ntags;++i)
146      {
147        m_offsetLookup.insert(DataMapType::value_type(tags[i],(i+1)*valsize));
148    }    }
149  }  }
150    
151    
152  DataTagged::DataTagged(const DataTagged& other)  DataTagged::DataTagged(const DataTagged& other)
153    : DataAbstract(other.getFunctionSpace()),    : parent(other.getFunctionSpace(),other.getShape()),
154    m_data(other.m_data),    m_offsetLookup(other.m_offsetLookup),
155    m_offsetLookup(other.m_offsetLookup)    m_data(other.m_data)
156  {  {
157    // copy constructor    // copy constructor
158    
159    // create the data view    // create the data view
160    DataArrayView temp(m_data,other.getPointDataView().getShape());  //   DataArrayView temp(m_data,other.getPointDataView().getShape());
161    setPointDataView(temp);  //   setPointDataView(temp);
162  }  }
163    
164  DataTagged::DataTagged(const DataConstant& other)  DataTagged::DataTagged(const DataConstant& other)
165    : DataAbstract(other.getFunctionSpace())    : parent(other.getFunctionSpace(),other.getShape())
166  {  {
167    // copy constructor    // copy constructor
168    
169      if (!other.getFunctionSpace().canTag())
170      {
171        throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
172      }
173    
174    // fill the default value with the constant value item from "other"    // fill the default value with the constant value item from "other"
175    const DataArrayView& value=other.getPointDataView();    int len = other.getNoValues();
   int len = value.noValues();  
176    m_data.resize(len,0.,len);    m_data.resize(len,0.,len);
177    for (int i=0; i<value.noValues(); i++) {    for (int i=0; i<len; i++) {
178      m_data[i]=value.getData(i);      m_data[i]=other.getVectorRO()[i];
179    }    }
180    }
181    
182    // create the data view  
183    DataArrayView temp(m_data,value.getShape());  // Create a new object by copying tags
184    setPointDataView(temp);  DataTagged::DataTagged(const FunctionSpace& what,
185                 const DataTypes::ShapeType& shape,
186             const DataTypes::ValueType& defaultvalue,
187                 const DataTagged* tagsource)
188     : parent(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*  DataAbstract*
223  DataTagged::getSlice(const DataArrayView::RegionType& region) const  DataTagged::deepCopy()
224    {
225      return new DataTagged(*this);
226    }
227    
228    DataAbstract*
229    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())    : parent(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    
   // create the data view  
   DataArrayView temp(m_data,regionShape);  
   setPointDataView(temp);  
   
249    // copy the default value from other to this    // copy the default value from other to this
250    getDefaultValue().copySlice(other.getDefaultValue(), regionLoopRange);    const DataTypes::ShapeType& otherShape=other.getShape();
251      const DataTypes::ValueType& otherData=other.getVectorRO();
252      DataTypes::copySlice(getVectorRW(),getShape(),getDefaultOffset(),otherData,otherShape,other.getDefaultOffset(), regionLoopRange);
253    
254    // loop through the tag values copying these    // loop through the tag values copying these
255    DataMapType::const_iterator pos;    DataMapType::const_iterator pos;
256    DataArrayView::ValueType::size_type tagOffset=getPointDataView().noValues();    DataTypes::ValueType::size_type tagOffset=getNoValues();
257    for (pos=other.m_offsetLookup.begin();pos!=other.m_offsetLookup.end();pos++){    for (pos=other.m_offsetLookup.begin();pos!=other.m_offsetLookup.end();pos++){
258      getPointDataView().copySlice(tagOffset,other.getPointDataView(),pos->second,regionLoopRange);      DataTypes::copySlice(m_data,getShape(),tagOffset,otherData, otherShape, pos->second, regionLoopRange);
259      m_offsetLookup.insert(DataMapType::value_type(pos->first,tagOffset));      m_offsetLookup.insert(DataMapType::value_type(pos->first,tagOffset));
260      tagOffset+=getPointDataView().noValues();      tagOffset+=getNoValues();
261    }    }
262  }  }
263    
264  void  void
265  DataTagged::setSlice(const DataAbstract* other,  DataTagged::setSlice(const DataAbstract* other,
266                       const DataArrayView::RegionType& region)                       const DataTypes::RegionType& region)
267  {  {
268    
269    // other must be another DataTagged object    // other must be another DataTagged object
# Line 184  DataTagged::setSlice(const DataAbstract* Line 273  DataTagged::setSlice(const DataAbstract*
273      throw DataException("Programming error - casting to DataTagged.");      throw DataException("Programming error - casting to DataTagged.");
274    }    }
275    
276      CHECK_FOR_EX_WRITE
277    
278    // determine shape of the specified region    // determine shape of the specified region
279    DataArrayView::ShapeType regionShape(DataArrayView::getResultSliceShape(region));    DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
280    
281    // modify region specification as needed to match rank of this object    // modify region specification as needed to match rank of this object
282    DataArrayView::RegionLoopRangeType regionLoopRange=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
283    
284    // ensure rank/shape of this object is compatible with specified region    // ensure rank/shape of this object is compatible with specified region
285    if (getPointDataView().getRank()!=region.size()) {    if (getRank()!=region.size()) {
286      throw DataException("Error - Invalid slice region.");      throw DataException("Error - Invalid slice region.");
287    }    }
288    if (otherTemp->getPointDataView().getRank()>0 && !other->getPointDataView().checkShape(regionShape)) {    if (otherTemp->getRank()>0 && !DataTypes::checkShape(other->getShape(),regionShape)) {
289      throw DataException (other->getPointDataView().createShapeErrorMessage(      throw DataException (DataTypes::createShapeErrorMessage(
290                           "Error - Couldn't copy slice due to shape mismatch.",regionShape));                           "Error - Couldn't copy slice due to shape mismatch.",regionShape,other->getShape()));
291    }    }
292    
293      const DataTypes::ValueType& otherData=otherTemp->getVectorRO();
294      const DataTypes::ShapeType& otherShape=otherTemp->getShape();
295    // copy slice from other default value to this default value    // copy slice from other default value to this default value
296    getDefaultValue().copySliceFrom(otherTemp->getDefaultValue(), regionLoopRange);    DataTypes::copySliceFrom(m_data,getShape(),getDefaultOffset(),otherData,otherShape,otherTemp->getDefaultOffset(),regionLoopRange);
297    
298    // 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
299    DataMapType::const_iterator pos;    DataMapType::const_iterator pos;
300    for (pos=otherTemp->m_offsetLookup.begin();pos!=otherTemp->m_offsetLookup.end();pos++) {    for (pos=otherTemp->m_offsetLookup.begin();pos!=otherTemp->m_offsetLookup.end();pos++) {
301      if (!isCurrentTag(pos->first)) {      if (!isCurrentTag(pos->first)) {
302        addTaggedValue(pos->first,getDefaultValue());        addTag(pos->first);
303      }      }
304    }    }
305    
306    // loop through the tag values copying slices from other to this    // loop through the tag values copying slices from other to this
307    for (pos=m_offsetLookup.begin();pos!=m_offsetLookup.end();pos++) {    for (pos=m_offsetLookup.begin();pos!=m_offsetLookup.end();pos++) {
308      getDataPointByTag(pos->first).copySliceFrom(otherTemp->getDataPointByTag(pos->first), regionLoopRange);      DataTypes::copySliceFrom(m_data,getShape(),getOffsetForTag(pos->first),otherData, otherShape, otherTemp->getOffsetForTag(pos->first), regionLoopRange);
309    
310    }    }
311    
312  }  }
# Line 248  DataTagged::getTagNumber(int dpno) Line 342  DataTagged::getTagNumber(int dpno)
342  }  }
343    
344  void  void
 DataTagged::setTaggedValues(const TagListType& tagKeys,  
                             const ValueListType& values)  
 {  
   addTaggedValues(tagKeys,values);  
 }  
   
 void  
345  DataTagged::setTaggedValue(int tagKey,  DataTagged::setTaggedValue(int tagKey,
346                             const DataArrayView& value)                 const DataTypes::ShapeType& pointshape,
347  {                             const ValueType& value,
348    if (!getPointDataView().checkShape(value.getShape())) {                 int dataOffset)
349        throw DataException(getPointDataView().createShapeErrorMessage(  {
350                            "Error - Cannot setTaggedValue due to shape mismatch.", value.getShape()));    if (!DataTypes::checkShape(getShape(), pointshape)) {
351          throw DataException(DataTypes::createShapeErrorMessage(
352                              "Error - Cannot setTaggedValue due to shape mismatch.", pointshape,getShape()));
353    }    }
354      CHECK_FOR_EX_WRITE
355    DataMapType::iterator pos(m_offsetLookup.find(tagKey));    DataMapType::iterator pos(m_offsetLookup.find(tagKey));
356    if (pos==m_offsetLookup.end()) {    if (pos==m_offsetLookup.end()) {
357      // tag couldn't be found so use addTaggedValue      // tag couldn't be found so use addTaggedValue
358      addTaggedValue(tagKey,value);      addTaggedValue(tagKey,pointshape, value, dataOffset);
359    } else {    } else {
360      // copy the values into the data array at the offset determined by m_offsetLookup      // copy the values into the data array at the offset determined by m_offsetLookup
361      int offset=pos->second;      int offset=pos->second;
362      for (int i=0; i<getPointDataView().noValues(); i++) {      for (unsigned int i=0; i<getNoValues(); i++) {
363        m_data[offset+i]=value.getData(i);        m_data[offset+i]=value[i+dataOffset];
364      }      }
365    }    }
366  }  }
367    
368    
369  void  void
370  DataTagged::addTaggedValues(const TagListType& tagKeys,  DataTagged::addTaggedValues(const TagListType& tagKeys,
371                              const ValueListType& values)                              const ValueBatchType& values,
372                                const ShapeType& vShape)
373  {  {
374      DataTypes::ValueType t(values.size(),0);
375      for (size_t i=0;i<values.size();++i)
376      {
377        t[i]=values[i];
378      }
379      addTaggedValues(tagKeys,t,vShape);
380    }
381    
382    
383    // Note: The check to see if vShape==our shape is done in the addTaggedValue method
384    void
385    DataTagged::addTaggedValues(const TagListType& tagKeys,
386                                const ValueType& values,
387                                const ShapeType& vShape)
388    {
389      unsigned int n=getNoValues();
390      unsigned int numVals=values.size()/n;
391    if (values.size()==0) {    if (values.size()==0) {
392      // copy the current default value for each of the tags      // copy the current default value for each of the tags
393      TagListType::const_iterator iT;      TagListType::const_iterator iT;
394      for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {      for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
395        // the point data view for DataTagged points at the default value        // the point data view for DataTagged points at the default value
396        addTaggedValue(*iT,getPointDataView());        addTag(*iT);
397      }      }
398    } else if (values.size()==1 && tagKeys.size()>1) {    } else if (numVals==1 && tagKeys.size()>1) {
399      // assume the one given value will be used for all tag values      // assume the one given value will be used for all tag values
400      TagListType::const_iterator iT;      TagListType::const_iterator iT;
401      for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {      for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
402        addTaggedValue(*iT,values[0]);        addTaggedValue(*iT, vShape, values,0);
403      }      }
404    } else {    } else {
405      if (tagKeys.size()!=values.size()) {      if (tagKeys.size()!=numVals) {
406        stringstream temp;        stringstream temp;
407        temp << "Error - (addTaggedValue) Number of tags: " << tagKeys.size()        temp << "Error - (addTaggedValue) Number of tags: " << tagKeys.size()
408         << " doesn't match number of values: " << values.size();         << " doesn't match number of values: " << values.size();
409        throw DataException(temp.str());        throw DataException(temp.str());
410      } else {      } else {
411        for (int i=0;i<tagKeys.size();i++) {        unsigned int i;
412          addTaggedValue(tagKeys[i],values[i]);        int offset=0;
413          for (i=0;i<tagKeys.size();i++ ,offset+=n) {
414            addTaggedValue(tagKeys[i],vShape,values,offset);
415        }        }
416      }      }
417    }    }
418  }  }
419    
420    
421    
422    
423  void  void
424  DataTagged::addTaggedValue(int tagKey,  DataTagged::addTaggedValue(int tagKey,
425                             const DataArrayView& value)                 const DataTypes::ShapeType& pointshape,
426  {                             const ValueType& value,
427    if (!getPointDataView().checkShape(value.getShape())) {                 int dataOffset)
428      throw DataException(getPointDataView().createShapeErrorMessage(  {
429                          "Error - Cannot addTaggedValue due to shape mismatch.", value.getShape()));    if (!DataTypes::checkShape(getShape(), pointshape)) {
430        throw DataException(DataTypes::createShapeErrorMessage(
431                            "Error - Cannot addTaggedValue due to shape mismatch.", pointshape,getShape()));
432    }    }
433      CHECK_FOR_EX_WRITE
434    DataMapType::iterator pos(m_offsetLookup.find(tagKey));    DataMapType::iterator pos(m_offsetLookup.find(tagKey));
435    if (pos!=m_offsetLookup.end()) {    if (pos!=m_offsetLookup.end()) {
436      // tag already exists so use setTaggedValue      // tag already exists so use setTaggedValue
437      setTaggedValue(tagKey,value);      setTaggedValue(tagKey,pointshape, value, dataOffset);
438    } else {    } else {
439      // 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
440      m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));      m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));
# Line 326  DataTagged::addTaggedValue(int tagKey, Line 443  DataTagged::addTaggedValue(int tagKey,
443      // 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
444      ValueType m_data_temp(m_data);      ValueType m_data_temp(m_data);
445      int oldSize=m_data.size();      int oldSize=m_data.size();
446      int newSize=m_data.size()+value.noValues();      int newSize=m_data.size()+getNoValues();
447      m_data.resize(newSize,0.,newSize);      m_data.resize(newSize,0.,newSize);
448      for (int i=0;i<oldSize;i++) {      for (int i=0;i<oldSize;i++) {
449        m_data[i]=m_data_temp[i];        m_data[i]=m_data_temp[i];
450      }      }
451      for (int i=0;i<value.noValues();i++) {      for (unsigned int i=0;i<getNoValues();i++) {
452        m_data[oldSize+i]=value.getData(i);        m_data[oldSize+i]=value[i+dataOffset];
453      }      }
454    }    }
455  }  }
456    
457    void
458    DataTagged::addTag(int tagKey)
459    {
460      CHECK_FOR_EX_WRITE
461      DataMapType::iterator pos(m_offsetLookup.find(tagKey));
462      if (pos!=m_offsetLookup.end()) {
463        // tag already exists so use setTaggedValue
464    //    setTaggedValue(tagKey,value);
465      } else {
466        // save the key and the location of its data in the lookup tab
467        m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));
468        // add the data given in "value" at the end of m_data
469        // need to make a temp copy of m_data, resize m_data, then copy
470        // all the old values plus the value to be added back into m_data
471        ValueType m_data_temp(m_data);
472        int oldSize=m_data.size();
473        int newSize=m_data.size()+getNoValues();
474        m_data.resize(newSize,0.,newSize);
475        for (int i=0;i<oldSize;i++) {
476          m_data[i]=m_data_temp[i];
477        }
478        for (unsigned int i=0;i<getNoValues();i++) {
479          m_data[oldSize+i]=m_data[m_defaultValueOffset+i];
480        }
481      }
482    }
483    
484    
485  double*  double*
486  DataTagged::getSampleDataByTag(int tag)  DataTagged::getSampleDataByTag(int tag)
487  {  {
488      CHECK_FOR_EX_WRITE
489    DataMapType::iterator pos(m_offsetLookup.find(tag));    DataMapType::iterator pos(m_offsetLookup.find(tag));
490    if (pos==m_offsetLookup.end()) {    if (pos==m_offsetLookup.end()) {
491      // tag couldn't be found so return the default value      // tag couldn't be found so return the default value
# Line 353  DataTagged::getSampleDataByTag(int tag) Line 499  DataTagged::getSampleDataByTag(int tag)
499  string  string
500  DataTagged::toString() const  DataTagged::toString() const
501  {  {
502      using namespace escript::DataTypes;
503      string empty="";
504    stringstream temp;    stringstream temp;
505    DataMapType::const_iterator i;    DataMapType::const_iterator i;
506    temp << "Tag(Default)" << endl;    temp << "Tag(Default)" << endl;
507    temp << getDefaultValue().toString() << endl;    temp << pointToString(m_data,getShape(),getDefaultOffset(),empty) << endl;
508    // create a temporary view as the offset will be changed    // create a temporary view as the offset will be changed
509    DataArrayView tempView(getPointDataView().getData(), getPointDataView().getShape());  //   DataArrayView tempView(getPointDataView().getData(), getPointDataView().getShape());
510    for (i=m_offsetLookup.begin();i!=m_offsetLookup.end();++i) {    for (i=m_offsetLookup.begin();i!=m_offsetLookup.end();++i) {
511      temp << "Tag(" << i->first << ")" << endl;      temp << "Tag(" << i->first << ")" << endl;
512      tempView.setOffset(i->second);      temp << pointToString(m_data,getShape(),i->second,empty) << endl;
513      temp << tempView.toString() << endl;  //     tempView.setOffset(i->second);
514    //     temp << tempView.toString() << endl;
515    }    }
516    return temp.str();    return temp.str();
517  }  }
518    
519  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
520  DataTagged::getPointOffset(int sampleNo,  DataTagged::getPointOffset(int sampleNo,
521                             int dataPointNo) const                             int dataPointNo) const
522  {  {
523    int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);    int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
524    DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));    DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));
525    DataArrayView::ValueType::size_type offset=m_defaultValueOffset;    DataTypes::ValueType::size_type offset=m_defaultValueOffset;
526    if (pos!=m_offsetLookup.end()) {    if (pos!=m_offsetLookup.end()) {
527      offset=pos->second;      offset=pos->second;
528    }    }
529    return offset;    return offset;
530  }  }
531    
532  DataArrayView  DataTypes::ValueType::size_type
533  DataTagged::getDataPointByTag(int tag) const  DataTagged::getPointOffset(int sampleNo,
534                               int dataPointNo)
535  {  {
536    DataMapType::const_iterator pos(m_offsetLookup.find(tag));    int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
537    DataArrayView::ValueType::size_type offset=m_defaultValueOffset;    DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));
538      DataTypes::ValueType::size_type offset=m_defaultValueOffset;
539    if (pos!=m_offsetLookup.end()) {    if (pos!=m_offsetLookup.end()) {
540      offset=pos->second;      offset=pos->second;
541    }    }
542    DataArrayView temp(getPointDataView());    return offset;
   temp.setOffset(offset);  
   return temp;  
543  }  }
544    
545  DataArrayView  DataTypes::ValueType::size_type
546  DataTagged::getDataPoint(int sampleNo,  DataTagged::getOffsetForTag(int tag) const
                          int dataPointNo)  
547  {  {
548    EsysAssert(validSampleNo(sampleNo),"(getDataPoint) Invalid sampleNo: " << sampleNo);    DataMapType::const_iterator pos(m_offsetLookup.find(tag));
549    int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);    DataTypes::ValueType::size_type offset=m_defaultValueOffset;
550    return getDataPointByTag(tagKey);    if (pos!=m_offsetLookup.end()) {
551        offset=pos->second;
552      }
553      return offset;
554  }  }
555    
556  int  DataTypes::ValueType::const_reference
557  DataTagged::archiveData(ofstream& archiveFile,  DataTagged::getDataByTagRO(int tag, DataTypes::ValueType::size_type i) const
                         const DataArrayView::ValueType::size_type noValues) const  
558  {  {
559    return(m_data.archiveData(archiveFile, noValues));    DataMapType::const_iterator pos(m_offsetLookup.find(tag));
560      DataTypes::ValueType::size_type offset=m_defaultValueOffset;
561      if (pos!=m_offsetLookup.end()) {
562        offset=pos->second;
563      }
564      return m_data[offset+i];
565  }  }
566    
567  int  DataTypes::ValueType::reference
568  DataTagged::extractData(ifstream& archiveFile,  DataTagged::getDataByTagRW(int tag, DataTypes::ValueType::size_type i)
                         const DataArrayView::ValueType::size_type noValues)  
569  {  {
570    return(m_data.extractData(archiveFile, noValues));    CHECK_FOR_EX_WRITE
571      DataMapType::const_iterator pos(m_offsetLookup.find(tag));
572      DataTypes::ValueType::size_type offset=m_defaultValueOffset;
573      if (pos!=m_offsetLookup.end()) {
574        offset=pos->second;
575      }
576      return m_data[offset+i];
577  }  }
578    
579  void  void
580  DataTagged::symmetric(DataAbstract* ev)  DataTagged::symmetric(DataAbstract* ev)
581  {  {
# Line 425  DataTagged::symmetric(DataAbstract* ev) Line 586  DataTagged::symmetric(DataAbstract* ev)
586    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
587    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
588    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
589      ValueType& evVec=temp_ev->getVectorRW();
590      const ShapeType& evShape=temp_ev->getShape();
591    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
592        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
593        DataArrayView thisView=getDataPointByTag(i->first);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
594        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
595        DataArrayView::symmetric(thisView,0,evView,0);        DataMaths::symmetric(m_data,getShape(),offset,evVec, evShape, evoffset);
596    }    }
597    DataArrayView::symmetric(getDefaultValue(),0,temp_ev->getDefaultValue(),0);    DataMaths::symmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
598  }  }
599    
600    
601  void  void
602  DataTagged::nonsymmetric(DataAbstract* ev)  DataTagged::nonsymmetric(DataAbstract* ev)
603  {  {
# Line 443  DataTagged::nonsymmetric(DataAbstract* e Line 608  DataTagged::nonsymmetric(DataAbstract* e
608    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
609    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
610    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
611      ValueType& evVec=temp_ev->getVectorRW();
612      const ShapeType& evShape=temp_ev->getShape();
613    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
614        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
615        DataArrayView thisView=getDataPointByTag(i->first);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
616        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
617        DataArrayView::nonsymmetric(thisView,0,evView,0);        DataMaths::nonsymmetric(m_data,getShape(),offset,evVec, evShape, evoffset);
618    }    }
619    DataArrayView::nonsymmetric(getDefaultValue(),0,temp_ev->getDefaultValue(),0);    DataMaths::nonsymmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
620  }  }
621    
622    
623  void  void
624  DataTagged::trace(DataAbstract* ev, int axis_offset)  DataTagged::trace(DataAbstract* ev, int axis_offset)
625  {  {
# Line 461  DataTagged::trace(DataAbstract* ev, int Line 630  DataTagged::trace(DataAbstract* ev, int
630    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
631    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
632    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
633      ValueType& evVec=temp_ev->getVectorRW();
634      const ShapeType& evShape=temp_ev->getShape();
635    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
636        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
637        DataArrayView thisView=getDataPointByTag(i->first);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
638        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
639        DataArrayView::trace(thisView,0,evView,0, axis_offset);        DataMaths::trace(m_data,getShape(),offset,evVec, evShape, evoffset, axis_offset);
640    }    }
641    DataArrayView::trace(getDefaultValue(),0,temp_ev->getDefaultValue(),0,axis_offset);    DataMaths::trace(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
642  }  }
643    
644  void  void
# Line 480  DataTagged::transpose(DataAbstract* ev, Line 651  DataTagged::transpose(DataAbstract* ev,
651    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
652    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
653    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
654      ValueType& evVec=temp_ev->getVectorRW();
655      const ShapeType& evShape=temp_ev->getShape();
656    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
657        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
658        DataArrayView thisView=getDataPointByTag(i->first);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
659        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
660        DataArrayView::transpose(thisView,0,evView,0, axis_offset);        DataMaths::transpose(m_data,getShape(),offset,evVec, evShape, evoffset, axis_offset);
661    }    }
662    DataArrayView::transpose(getDefaultValue(),0,temp_ev->getDefaultValue(),0,axis_offset);    DataMaths::transpose(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
663  }  }
664    
665  void  void
# Line 499  DataTagged::swapaxes(DataAbstract* ev, i Line 672  DataTagged::swapaxes(DataAbstract* ev, i
672    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
673    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
674    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
675      ValueType& evVec=temp_ev->getVectorRW();
676      const ShapeType& evShape=temp_ev->getShape();
677    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
678        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
679        DataArrayView thisView=getDataPointByTag(i->first);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
680        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
681        DataArrayView::swapaxes(thisView,0,evView,0,axis0,axis1);        DataMaths::swapaxes(m_data,getShape(),offset,evVec, evShape, evoffset,axis0,axis1);
682    }    }
683    DataArrayView::swapaxes(getDefaultValue(),0,temp_ev->getDefaultValue(),0,axis0,axis1);    DataMaths::swapaxes(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis0,axis1);
684  }  }
685    
686  void  void
# Line 518  DataTagged::eigenvalues(DataAbstract* ev Line 693  DataTagged::eigenvalues(DataAbstract* ev
693    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
694    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
695    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
696      ValueType& evVec=temp_ev->getVectorRW();
697      const ShapeType& evShape=temp_ev->getShape();
698    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
699        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
700        DataArrayView thisView=getDataPointByTag(i->first);  //       DataArrayView thisView=getDataPointByTag(i->first);
701        DataArrayView evView=temp_ev->getDataPointByTag(i->first);  //       DataArrayView evView=temp_ev->getDataPointByTag(i->first);
702        DataArrayView::eigenvalues(thisView,0,evView,0);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
703          DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
704          DataMaths::eigenvalues(m_data,getShape(),offset,evVec, evShape, evoffset);
705    }    }
706    DataArrayView::eigenvalues(getDefaultValue(),0,temp_ev->getDefaultValue(),0);    DataMaths::eigenvalues(m_data,getShape(),getDefaultOffset(),evVec, evShape, temp_ev->getDefaultOffset());
707  }  }
708  void  void
709  DataTagged::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)  DataTagged::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
# Line 540  DataTagged::eigenvalues_and_eigenvectors Line 719  DataTagged::eigenvalues_and_eigenvectors
719    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
720    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
721    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
722      ValueType& evVec=temp_ev->getVectorRW();
723      const ShapeType& evShape=temp_ev->getShape();
724      ValueType& VVec=temp_V->getVectorRW();
725      const ShapeType& VShape=temp_V->getShape();
726    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
727        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
728        temp_V->addTaggedValue(i->first,temp_V->getDefaultValue());        temp_V->addTag(i->first);
729        DataArrayView thisView=getDataPointByTag(i->first);  /*      DataArrayView thisView=getDataPointByTag(i->first);
730        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataArrayView evView=temp_ev->getDataPointByTag(i->first);
731        DataArrayView VView=temp_V->getDataPointByTag(i->first);        DataArrayView VView=temp_V->getDataPointByTag(i->first);*/
732        DataArrayView::eigenvalues_and_eigenvectors(thisView,0,evView,0,VView,0,tol);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
733          DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
734          DataTypes::ValueType::size_type Voffset=temp_V->getOffsetForTag(i->first);
735    /*      DataArrayView::eigenvalues_and_eigenvectors(thisView,0,evView,0,VView,0,tol);*/
736          DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),offset,evVec, evShape, evoffset,VVec,VShape,Voffset,tol);
737    
738    }    }
739    DataArrayView::eigenvalues_and_eigenvectors(getDefaultValue(),0,    DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),getDefaultOffset(),evVec, evShape,
740                                                temp_ev->getDefaultValue(),0,                        temp_ev->getDefaultOffset(),VVec,VShape,
741                                                temp_V->getDefaultValue(),0,                        temp_V->getDefaultOffset(), tol);
742                                                tol);  
743    
744    }
745    
746    void
747    DataTagged::setToZero(){
748        CHECK_FOR_EX_WRITE
749        DataTypes::ValueType::size_type n=m_data.size();
750        for (int i=0; i<n ;++i) m_data[i]=0.;
751  }  }
752    
753  void  void
754  DataTagged::dump(const std::string fileName) const  DataTagged::dump(const std::string fileName) const
755  {  {
    #ifdef PASO_MPI  
    throw DataException("Error - DataTagged:: dump is not implemented for MPI yet.")  
    #endif  
756     #ifdef USE_NETCDF     #ifdef USE_NETCDF
757     const int ldims=DataArrayView::maxRank+1;     const int ldims=DataTypes::maxRank+1;
758     const NcDim* ncdims[ldims];     const NcDim* ncdims[ldims];
759     NcVar *var, *tags_var;     NcVar *var, *tags_var;
760     int rank = getPointDataView().getRank();     int rank = getRank();
761     int type=  getFunctionSpace().getTypeCode();     int type=  getFunctionSpace().getTypeCode();
762     int ndims =0;     int ndims =0;
763     long dims[ldims];     long dims[ldims];
764     DataArrayView::ShapeType shape = getPointDataView().getShape();     const double* d_ptr=&(m_data[0]);
765       DataTypes::ShapeType shape = getShape();
766       int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
767       int mpi_num=getFunctionSpace().getDomain()->getMPISize();
768    #ifdef PASO_MPI
769       MPI_Status status;
770    #endif
771    
772    #ifdef PASO_MPI
773       /* Serialize NetCDF I/O */
774       if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81803, MPI_COMM_WORLD, &status);
775    #endif
776    
777     // netCDF error handler     // netCDF error handler
778     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
779     // Create the file.     // Create the file.
780     NcFile dataFile(fileName.c_str(), NcFile::Replace);     char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
781       NcFile dataFile(newFileName, NcFile::Replace);
782     // check if writing was successful     // check if writing was successful
783     if (!dataFile.is_valid())     if (!dataFile.is_valid())
784          throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");          throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");
785     if (!dataFile.add_att("type","tagged") )     if (!dataFile.add_att("type_id",1) )
786          throw DataException("Error - DataTagged:: appending data type to netCDF file failed.");          throw DataException("Error - DataTagged:: appending data type to netCDF file failed.");
787     if (!dataFile.add_att("rank",rank) )     if (!dataFile.add_att("rank",rank) )
788          throw DataException("Error - DataTagged:: appending rank attribute to netCDF file failed.");          throw DataException("Error - DataTagged:: appending rank attribute to netCDF file failed.");
# Line 589  DataTagged::dump(const std::string fileN Line 792  DataTagged::dump(const std::string fileN
792     if ( rank >0 ) {     if ( rank >0 ) {
793         dims[0]=shape[0];         dims[0]=shape[0];
794         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
795              throw DataException("Error - DataTagged:: appending ncdimsion 0 to netCDF file failed.");              throw DataException("Error - DataTagged:: appending ncdimension 0 to netCDF file failed.");
796     }     }
797     if ( rank >1 ) {     if ( rank >1 ) {
798         dims[1]=shape[1];         dims[1]=shape[1];
799         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
800              throw DataException("Error - DataTagged:: appending ncdimsion 1 to netCDF file failed.");              throw DataException("Error - DataTagged:: appending ncdimension 1 to netCDF file failed.");
801     }     }
802     if ( rank >2 ) {     if ( rank >2 ) {
803         dims[2]=shape[2];         dims[2]=shape[2];
804         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
805              throw DataException("Error - DataTagged:: appending ncdimsion 2 to netCDF file failed.");              throw DataException("Error - DataTagged:: appending ncdimension 2 to netCDF file failed.");
806     }     }
807     if ( rank >3 ) {     if ( rank >3 ) {
808         dims[3]=shape[3];         dims[3]=shape[3];
809         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
810              throw DataException("Error - DataTagged:: appending ncdimsion 3 to netCDF file failed.");              throw DataException("Error - DataTagged:: appending ncdimension 3 to netCDF file failed.");
811     }     }
812     const DataTagged::DataMapType& thisLookup=getTagLookup();     const DataTagged::DataMapType& thisLookup=getTagLookup();
813     DataTagged::DataMapType::const_iterator i;     DataTagged::DataMapType::const_iterator i;
814     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
815     int ntags=1;     int ntags=1;
816     for (i=thisLookup.begin();i!=thisLookupEnd;i++) ntags++;     for (i=thisLookup.begin();i!=thisLookupEnd;i++) ntags++;
817     int* tags =(int*) malloc(ntags*sizeof(int));     int* tags =(int*) esysUtils::malloc(ntags*sizeof(int));
818     int c=1;     int c=1;
819     tags[0]=-1;     tags[0]=-1;
820     for (i=thisLookup.begin();i!=thisLookupEnd;i++) tags[c++]=i->first;     for (i=thisLookup.begin();i!=thisLookupEnd;i++) tags[c++]=i->first;
821     dims[rank]=ntags;     dims[rank]=ntags;
822     if (! (ncdims[rank] = dataFile.add_dim("num_tags", dims[rank])) )     if (! (ncdims[rank] = dataFile.add_dim("num_tags", dims[rank])) )
823     {     {
824         free(tags);         esysUtils::free(tags);
825             throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");             throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");
826     }     }
827     if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )     if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )
828     {     {
829      free(tags);      esysUtils::free(tags);
830          throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");          throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");
831     }     }
832     if (! (tags_var->put(tags,dims[rank])) )     if (! (tags_var->put(tags,dims[rank])) )
833     {     {
834      free(tags);      esysUtils::free(tags);
835          throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");          throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");
836     }     }
837     if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )     if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
838     {     {
839      free(tags);      esysUtils::free(tags);
840          throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");          throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");
841     }     }
842     if (! (var->put(&m_data[0],dims)) )     if (! (var->put(d_ptr,dims)) )
843     {     {
844      free(tags);      esysUtils::free(tags);
845          throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");          throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");
846     }     }
847    #ifdef PASO_MPI
848       if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81803, MPI_COMM_WORLD);
849    #endif
850     #else     #else
851     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.");
852     #endif     #endif
853  }  }
854    
855    DataTypes::ValueType&
856    DataTagged::getVectorRW()
857    {
858        CHECK_FOR_EX_WRITE
859        return m_data;
860    }
861    
862    const DataTypes::ValueType&
863    DataTagged::getVectorRO() const
864    {
865        return m_data;
866    }
867    
868  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1023  
changed lines
  Added in v.2548

  ViewVC Help
Powered by ViewVC 1.1.26