/[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 1455 by phornby, Thu Feb 28 17:19:44 2008 UTC revision 2766 by jfenwick, Mon Nov 23 05:10:53 2009 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2009 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    
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);
   DataArrayView temp(m_data,DataArrayView::ShapeType());  
   setPointDataView(temp);  
 }  
   
 DataTagged::DataTagged(const TagListType& tagKeys,  
                const ValueListType& values,  
                const DataArrayView& defaultValue,  
                const FunctionSpace& what)  
   : DataAbstract(what)  
 {  
   // constructor  
   
   // initialise the array of data values  
   // the default value is always the first item in the values list  
   int len = defaultValue.noValues();  
   m_data.resize(len,0.,len);  
   for (int i=0; i<defaultValue.noValues(); i++) {  
     m_data[i]=defaultValue.getData(i);  
   }  
   
   // create the data view  
   DataArrayView temp(m_data,defaultValue.getShape());  
   setPointDataView(temp);  
   
   // add remaining tags and values  
   addTaggedValues(tagKeys,values);  
45  }  }
46    
47  DataTagged::DataTagged(const FunctionSpace& what,  DataTagged::DataTagged(const FunctionSpace& what,
48                         const DataArrayView::ShapeType &shape,                         const DataTypes::ShapeType &shape,
49                         const int tags[],                         const int tags[],
50                         const ValueType& data)                         const ValueType& data)
51    : DataAbstract(what)    : parent(what,shape)
52  {  {
53    // alternative constructor    // alternative constructor
54    // not unit_tested tested yet    // not unit_tested tested yet
55      // It is not explicitly unit tested yet, but it is called from DataFactory
56    
57      if (!what.canTag())
58      {
59        throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
60      }
61    // copy the data    // copy the data
62    m_data=data;    m_data=data;
63    
64    // create the view of the data    // we can't rely on the tag array to give us the number of tags so
65    DataArrayView tempView(m_data,shape);    // use the data we have been passed
66    setPointDataView(tempView);    int valsize=DataTypes::noValues(shape);
67      int ntags=data.size()/valsize;
68    
69    // create the tag lookup map    // create the tag lookup map
70    for (int sampleNo=0; sampleNo<getNumSamples(); sampleNo++) {    // we assume that the first value and first tag are the default value so we skip
71      m_offsetLookup.insert(DataMapType::value_type(sampleNo,tags[sampleNo]));    for (int i=1;i<ntags;++i)
72      {
73        m_offsetLookup.insert(DataMapType::value_type(tags[i],i*valsize));
74    }    }
75  }  }
76    
77  DataTagged::DataTagged(const FunctionSpace& what,  DataTagged::DataTagged(const FunctionSpace& what,
78                         const DataArrayView::ShapeType &shape,                         const DataTypes::ShapeType &shape,
79                         const TagListType& tags,                         const TagListType& tags,
80                         const ValueType& data)                         const ValueType& data)
81    : DataAbstract(what)    : parent(what,shape)
82  {  {
83    // alternative constructor    // alternative constructor
84    // not unit_tested tested yet  
85      if (!what.canTag())
86      {
87        throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
88      }
89    
90    // copy the data    // copy the data
91    m_data=data;    m_data=data;
92    
93    // create the view of the data    // create the view of the data
94    DataArrayView tempView(m_data,shape);  //   DataArrayView tempView(m_data,shape);
95    setPointDataView(tempView);  //   setPointDataView(tempView);
96    
97      // create the tag lookup map
98    
99    //   for (int sampleNo=0; sampleNo<getNumSamples(); sampleNo++) {
100    //     m_offsetLookup.insert(DataMapType::value_type(sampleNo,tags[sampleNo]));
101    //   }
102    
103      // The above code looks like it will create a map the wrong way around
104    
105      int valsize=DataTypes::noValues(shape);
106      int npoints=(data.size()/valsize)-1;
107      int ntags=tags.size();
108      if (ntags>npoints)
109      {     // This throw is not unit tested yet
110        throw DataException("Programming error - Too many tags for the supplied values.");
111      }
112    
113    // create the tag lookup map    // create the tag lookup map
114    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)
115      m_offsetLookup.insert(DataMapType::value_type(sampleNo,tags[sampleNo]));    for (int i=0;i<ntags;++i)
116      {
117        m_offsetLookup.insert(DataMapType::value_type(tags[i],(i+1)*valsize));
118    }    }
119  }  }
120    
121    
122  DataTagged::DataTagged(const DataTagged& other)  DataTagged::DataTagged(const DataTagged& other)
123    : DataAbstract(other.getFunctionSpace()),    : parent(other.getFunctionSpace(),other.getShape()),
124    m_data(other.m_data),    m_offsetLookup(other.m_offsetLookup),
125    m_offsetLookup(other.m_offsetLookup)    m_data(other.m_data)
126  {  {
127    // copy constructor    // copy constructor
128    
129    // create the data view    // create the data view
130    DataArrayView temp(m_data,other.getPointDataView().getShape());  //   DataArrayView temp(m_data,other.getPointDataView().getShape());
131    setPointDataView(temp);  //   setPointDataView(temp);
132  }  }
133    
134  DataTagged::DataTagged(const DataConstant& other)  DataTagged::DataTagged(const DataConstant& other)
135    : DataAbstract(other.getFunctionSpace())    : parent(other.getFunctionSpace(),other.getShape())
136  {  {
137    // copy constructor    // copy constructor
138    
139      if (!other.getFunctionSpace().canTag())
140      {
141        throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
142      }
143    
144    // fill the default value with the constant value item from "other"    // fill the default value with the constant value item from "other"
145    const DataArrayView& value=other.getPointDataView();    int len = other.getNoValues();
   int len = value.noValues();  
146    m_data.resize(len,0.,len);    m_data.resize(len,0.,len);
147    for (int i=0; i<value.noValues(); i++) {    for (int i=0; i<len; i++) {
148      m_data[i]=value.getData(i);      m_data[i]=other.getVectorRO()[i];
149    }    }
150    }
151    
152    // create the data view  
153    DataArrayView temp(m_data,value.getShape());  // Create a new object by copying tags
154    setPointDataView(temp);  DataTagged::DataTagged(const FunctionSpace& what,
155                 const DataTypes::ShapeType& shape,
156             const DataTypes::ValueType& defaultvalue,
157                 const DataTagged* tagsource)
158     : parent(what,shape)
159    {
160    // This constructor has not been unit tested yet
161    
162      if (defaultvalue.size()!=DataTypes::noValues(shape)) {
163        throw DataException("Programming error - defaultvalue does not match supplied shape.");
164      }
165    
166    
167      if (!what.canTag())
168      {
169        throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
170      }
171    
172      if (tagsource!=0)
173      {
174           m_data.resize(defaultvalue.size(),0.);   // since this is tagged data, we should have blocksize=1
175    
176           DataTagged::DataMapType::const_iterator i;
177           for (i=tagsource->getTagLookup().begin();i!=tagsource->getTagLookup().end();i++) {
178          addTag(i->first);
179           }
180      }
181      else
182      {
183        m_data.resize(defaultvalue.size());
184      }
185    
186      // need to set the default value ....
187      for (int i=0; i<defaultvalue.size(); i++) {
188         m_data[i]=defaultvalue[i];
189      }
190    }
191    
192    DataAbstract*
193    DataTagged::deepCopy()
194    {
195      return new DataTagged(*this);
196  }  }
197    
198  DataAbstract*  DataAbstract*
199  DataTagged::getSlice(const DataArrayView::RegionType& region) const  DataTagged::getSlice(const DataTypes::RegionType& region) const
200  {  {
201    return new DataTagged(*this, region);    return new DataTagged(*this, region);
202  }  }
203    
204  DataTagged::DataTagged(const DataTagged& other,  DataTagged::DataTagged(const DataTagged& other,
205                 const DataArrayView::RegionType& region)                 const DataTypes::RegionType& region)
206    : DataAbstract(other.getFunctionSpace())    : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
207  {  {
208    // slice constructor    // slice constructor
209    
210    // get the shape of the slice to copy from other    // get the shape of the slice to copy from other
211    DataArrayView::ShapeType regionShape(DataArrayView::getResultSliceShape(region));    DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
212    DataArrayView::RegionLoopRangeType regionLoopRange=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
213    
214    // allocate enough space in this for all values    // allocate enough space in this for all values
215    // (need to add one to allow for the default value)    // (need to add one to allow for the default value)
216    int len = DataArrayView::noValues(regionShape)*(other.m_offsetLookup.size()+1);    int len = DataTypes::noValues(regionShape)*(other.m_offsetLookup.size()+1);
217    m_data.resize(len,0.0,len);    m_data.resize(len,0.0,len);
218    
   // create the data view  
   DataArrayView temp(m_data,regionShape);  
   setPointDataView(temp);  
   
219    // copy the default value from other to this    // copy the default value from other to this
220    getDefaultValue().copySlice(other.getDefaultValue(), regionLoopRange);    const DataTypes::ShapeType& otherShape=other.getShape();
221      const DataTypes::ValueType& otherData=other.getVectorRO();
222      DataTypes::copySlice(getVectorRW(),getShape(),getDefaultOffset(),otherData,otherShape,other.getDefaultOffset(), regionLoopRange);
223    
224    // loop through the tag values copying these    // loop through the tag values copying these
225    DataMapType::const_iterator pos;    DataMapType::const_iterator pos;
226    DataArrayView::ValueType::size_type tagOffset=getPointDataView().noValues();    DataTypes::ValueType::size_type tagOffset=getNoValues();
227    for (pos=other.m_offsetLookup.begin();pos!=other.m_offsetLookup.end();pos++){    for (pos=other.m_offsetLookup.begin();pos!=other.m_offsetLookup.end();pos++){
228      getPointDataView().copySlice(tagOffset,other.getPointDataView(),pos->second,regionLoopRange);      DataTypes::copySlice(m_data,getShape(),tagOffset,otherData, otherShape, pos->second, regionLoopRange);
229      m_offsetLookup.insert(DataMapType::value_type(pos->first,tagOffset));      m_offsetLookup.insert(DataMapType::value_type(pos->first,tagOffset));
230      tagOffset+=getPointDataView().noValues();      tagOffset+=getNoValues();
231    }    }
232  }  }
233    
234  void  void
235  DataTagged::setSlice(const DataAbstract* other,  DataTagged::setSlice(const DataAbstract* other,
236                       const DataArrayView::RegionType& region)                       const DataTypes::RegionType& region)
237  {  {
238    
239    // other must be another DataTagged object    // other must be another DataTagged object
# Line 185  DataTagged::setSlice(const DataAbstract* Line 243  DataTagged::setSlice(const DataAbstract*
243      throw DataException("Programming error - casting to DataTagged.");      throw DataException("Programming error - casting to DataTagged.");
244    }    }
245    
246      CHECK_FOR_EX_WRITE
247    
248    // determine shape of the specified region    // determine shape of the specified region
249    DataArrayView::ShapeType regionShape(DataArrayView::getResultSliceShape(region));    DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
250    
251    // modify region specification as needed to match rank of this object    // modify region specification as needed to match rank of this object
252    DataArrayView::RegionLoopRangeType regionLoopRange=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
253    
254    // ensure rank/shape of this object is compatible with specified region    // ensure rank/shape of this object is compatible with specified region
255    if (getPointDataView().getRank()!=region.size()) {    if (getRank()!=region.size()) {
256      throw DataException("Error - Invalid slice region.");      throw DataException("Error - Invalid slice region.");
257    }    }
258    if (otherTemp->getPointDataView().getRank()>0 && !other->getPointDataView().checkShape(regionShape)) {    if (otherTemp->getRank()>0 && !DataTypes::checkShape(other->getShape(),regionShape)) {
259      throw DataException (other->getPointDataView().createShapeErrorMessage(      throw DataException (DataTypes::createShapeErrorMessage(
260                           "Error - Couldn't copy slice due to shape mismatch.",regionShape));                           "Error - Couldn't copy slice due to shape mismatch.",regionShape,other->getShape()));
261    }    }
262    
263      const DataTypes::ValueType& otherData=otherTemp->getVectorRO();
264      const DataTypes::ShapeType& otherShape=otherTemp->getShape();
265    // copy slice from other default value to this default value    // copy slice from other default value to this default value
266    getDefaultValue().copySliceFrom(otherTemp->getDefaultValue(), regionLoopRange);    DataTypes::copySliceFrom(m_data,getShape(),getDefaultOffset(),otherData,otherShape,otherTemp->getDefaultOffset(),regionLoopRange);
267    
268    // 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
269    DataMapType::const_iterator pos;    DataMapType::const_iterator pos;
270    for (pos=otherTemp->m_offsetLookup.begin();pos!=otherTemp->m_offsetLookup.end();pos++) {    for (pos=otherTemp->m_offsetLookup.begin();pos!=otherTemp->m_offsetLookup.end();pos++) {
271      if (!isCurrentTag(pos->first)) {      if (!isCurrentTag(pos->first)) {
272        addTaggedValue(pos->first,getDefaultValue());        addTag(pos->first);
273      }      }
274    }    }
275    
276    // loop through the tag values copying slices from other to this    // loop through the tag values copying slices from other to this
277    for (pos=m_offsetLookup.begin();pos!=m_offsetLookup.end();pos++) {    for (pos=m_offsetLookup.begin();pos!=m_offsetLookup.end();pos++) {
278      getDataPointByTag(pos->first).copySliceFrom(otherTemp->getDataPointByTag(pos->first), regionLoopRange);      DataTypes::copySliceFrom(m_data,getShape(),getOffsetForTag(pos->first),otherData, otherShape, otherTemp->getOffsetForTag(pos->first), regionLoopRange);
279    
280    }    }
281    
282  }  }
# Line 249  DataTagged::getTagNumber(int dpno) Line 312  DataTagged::getTagNumber(int dpno)
312  }  }
313    
314  void  void
 DataTagged::setTaggedValues(const TagListType& tagKeys,  
                             const ValueListType& values)  
 {  
   addTaggedValues(tagKeys,values);  
 }  
   
 void  
315  DataTagged::setTaggedValue(int tagKey,  DataTagged::setTaggedValue(int tagKey,
316                             const DataArrayView& value)                 const DataTypes::ShapeType& pointshape,
317  {                             const ValueType& value,
318    if (!getPointDataView().checkShape(value.getShape())) {                 int dataOffset)
319        throw DataException(getPointDataView().createShapeErrorMessage(  {
320                            "Error - Cannot setTaggedValue due to shape mismatch.", value.getShape()));    if (!DataTypes::checkShape(getShape(), pointshape)) {
321          throw DataException(DataTypes::createShapeErrorMessage(
322                              "Error - Cannot setTaggedValue due to shape mismatch.", pointshape,getShape()));
323    }    }
324      CHECK_FOR_EX_WRITE
325    DataMapType::iterator pos(m_offsetLookup.find(tagKey));    DataMapType::iterator pos(m_offsetLookup.find(tagKey));
326    if (pos==m_offsetLookup.end()) {    if (pos==m_offsetLookup.end()) {
327      // tag couldn't be found so use addTaggedValue      // tag couldn't be found so use addTaggedValue
328      addTaggedValue(tagKey,value);      addTaggedValue(tagKey,pointshape, value, dataOffset);
329    } else {    } else {
330      // 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
331      int offset=pos->second;      int offset=pos->second;
332      for (int i=0; i<getPointDataView().noValues(); i++) {      for (unsigned int i=0; i<getNoValues(); i++) {
333        m_data[offset+i]=value.getData(i);        m_data[offset+i]=value[i+dataOffset];
334      }      }
335    }    }
336  }  }
337    
338    
339  void  void
340  DataTagged::addTaggedValues(const TagListType& tagKeys,  DataTagged::addTaggedValues(const TagListType& tagKeys,
341                              const ValueListType& values)                              const ValueBatchType& values,
342                                const ShapeType& vShape)
343  {  {
344      DataTypes::ValueType t(values.size(),0);
345      for (size_t i=0;i<values.size();++i)
346      {
347        t[i]=values[i];
348      }
349      addTaggedValues(tagKeys,t,vShape);
350    }
351    
352    
353    // Note: The check to see if vShape==our shape is done in the addTaggedValue method
354    void
355    DataTagged::addTaggedValues(const TagListType& tagKeys,
356                                const ValueType& values,
357                                const ShapeType& vShape)
358    {
359      unsigned int n=getNoValues();
360      unsigned int numVals=values.size()/n;
361    if (values.size()==0) {    if (values.size()==0) {
362      // copy the current default value for each of the tags      // copy the current default value for each of the tags
363      TagListType::const_iterator iT;      TagListType::const_iterator iT;
364      for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {      for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
365        // the point data view for DataTagged points at the default value        // the point data view for DataTagged points at the default value
366        addTaggedValue(*iT,getPointDataView());        addTag(*iT);
367      }      }
368    } else if (values.size()==1 && tagKeys.size()>1) {    } else if (numVals==1 && tagKeys.size()>1) {
369      // assume the one given value will be used for all tag values      // assume the one given value will be used for all tag values
370      TagListType::const_iterator iT;      TagListType::const_iterator iT;
371      for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {      for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
372        addTaggedValue(*iT,values[0]);        addTaggedValue(*iT, vShape, values,0);
373      }      }
374    } else {    } else {
375      if (tagKeys.size()!=values.size()) {      if (tagKeys.size()!=numVals) {
376        stringstream temp;        stringstream temp;
377        temp << "Error - (addTaggedValue) Number of tags: " << tagKeys.size()        temp << "Error - (addTaggedValue) Number of tags: " << tagKeys.size()
378         << " doesn't match number of values: " << values.size();         << " doesn't match number of values: " << values.size();
379        throw DataException(temp.str());        throw DataException(temp.str());
380      } else {      } else {
381        unsigned int i;        unsigned int i;
382        for (i=0;i<tagKeys.size();i++) {        int offset=0;
383          addTaggedValue(tagKeys[i],values[i]);        for (i=0;i<tagKeys.size();i++ ,offset+=n) {
384            addTaggedValue(tagKeys[i],vShape,values,offset);
385        }        }
386      }      }
387    }    }
388  }  }
389    
390    
391    
392    
393  void  void
394  DataTagged::addTaggedValue(int tagKey,  DataTagged::addTaggedValue(int tagKey,
395                             const DataArrayView& value)                 const DataTypes::ShapeType& pointshape,
396  {                             const ValueType& value,
397    if (!getPointDataView().checkShape(value.getShape())) {                 int dataOffset)
398      throw DataException(getPointDataView().createShapeErrorMessage(  {
399                          "Error - Cannot addTaggedValue due to shape mismatch.", value.getShape()));    if (!DataTypes::checkShape(getShape(), pointshape)) {
400        throw DataException(DataTypes::createShapeErrorMessage(
401                            "Error - Cannot addTaggedValue due to shape mismatch.", pointshape,getShape()));
402    }    }
403      CHECK_FOR_EX_WRITE
404    DataMapType::iterator pos(m_offsetLookup.find(tagKey));    DataMapType::iterator pos(m_offsetLookup.find(tagKey));
405    if (pos!=m_offsetLookup.end()) {    if (pos!=m_offsetLookup.end()) {
406      // tag already exists so use setTaggedValue      // tag already exists so use setTaggedValue
407      setTaggedValue(tagKey,value);      setTaggedValue(tagKey,pointshape, value, dataOffset);
408    } else {    } else {
409      // 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
410      m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));      m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));
# Line 328  DataTagged::addTaggedValue(int tagKey, Line 413  DataTagged::addTaggedValue(int tagKey,
413      // 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
414      ValueType m_data_temp(m_data);      ValueType m_data_temp(m_data);
415      int oldSize=m_data.size();      int oldSize=m_data.size();
416      int newSize=m_data.size()+value.noValues();      int newSize=m_data.size()+getNoValues();
417      m_data.resize(newSize,0.,newSize);      m_data.resize(newSize,0.,newSize);
418      for (int i=0;i<oldSize;i++) {      for (int i=0;i<oldSize;i++) {
419        m_data[i]=m_data_temp[i];        m_data[i]=m_data_temp[i];
420      }      }
421      for (int i=0;i<value.noValues();i++) {      for (unsigned int i=0;i<getNoValues();i++) {
422        m_data[oldSize+i]=value.getData(i);        m_data[oldSize+i]=value[i+dataOffset];
423      }      }
424    }    }
425  }  }
426    
427    void
428    DataTagged::addTag(int tagKey)
429    {
430      CHECK_FOR_EX_WRITE
431      DataMapType::iterator pos(m_offsetLookup.find(tagKey));
432      if (pos!=m_offsetLookup.end()) {
433        // tag already exists so use setTaggedValue
434    //    setTaggedValue(tagKey,value);
435      } else {
436        // save the key and the location of its data in the lookup tab
437        m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));
438        // add the data given in "value" at the end of m_data
439        // need to make a temp copy of m_data, resize m_data, then copy
440        // all the old values plus the value to be added back into m_data
441        ValueType m_data_temp(m_data);
442        int oldSize=m_data.size();
443        int newSize=m_data.size()+getNoValues();
444        m_data.resize(newSize,0.,newSize);
445        for (int i=0;i<oldSize;i++) {
446          m_data[i]=m_data_temp[i];
447        }
448        for (unsigned int i=0;i<getNoValues();i++) {
449          m_data[oldSize+i]=m_data[m_defaultValueOffset+i];
450        }
451      }
452    }
453    
454    
455  double*  double*
456  DataTagged::getSampleDataByTag(int tag)  DataTagged::getSampleDataByTag(int tag)
457  {  {
458      CHECK_FOR_EX_WRITE
459    DataMapType::iterator pos(m_offsetLookup.find(tag));    DataMapType::iterator pos(m_offsetLookup.find(tag));
460    if (pos==m_offsetLookup.end()) {    if (pos==m_offsetLookup.end()) {
461      // tag couldn't be found so return the default value      // tag couldn't be found so return the default value
# Line 352  DataTagged::getSampleDataByTag(int tag) Line 466  DataTagged::getSampleDataByTag(int tag)
466    }    }
467  }  }
468    
469    
470    bool
471    DataTagged::hasNaN() const
472    {
473        for (ValueType::size_type i=0;i<m_data.size();++i)
474        {
475            if (nancheck(m_data[i]))    // can't assume we have new standard NaN checking
476            {
477                return true;
478            }
479        }
480        return false;
481    }
482    
483  string  string
484  DataTagged::toString() const  DataTagged::toString() const
485  {  {
486      using namespace escript::DataTypes;
487      string empty="";
488    stringstream temp;    stringstream temp;
489    DataMapType::const_iterator i;    DataMapType::const_iterator i;
490    temp << "Tag(Default)" << endl;    temp << "Tag(Default)" << endl;
491    temp << getDefaultValue().toString() << endl;    temp << pointToString(m_data,getShape(),getDefaultOffset(),empty) << endl;
492    // create a temporary view as the offset will be changed    // create a temporary view as the offset will be changed
493    DataArrayView tempView(getPointDataView().getData(), getPointDataView().getShape());  //   DataArrayView tempView(getPointDataView().getData(), getPointDataView().getShape());
494    for (i=m_offsetLookup.begin();i!=m_offsetLookup.end();++i) {    for (i=m_offsetLookup.begin();i!=m_offsetLookup.end();++i) {
495      temp << "Tag(" << i->first << ")" << endl;      temp << "Tag(" << i->first << ")" << endl;
496      tempView.setOffset(i->second);      temp << pointToString(m_data,getShape(),i->second,empty) << endl;
497      temp << tempView.toString() << endl;  //     tempView.setOffset(i->second);
498    //     temp << tempView.toString() << endl;
499    }    }
500    return temp.str();    return temp.str();
501  }  }
502    
503  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
504  DataTagged::getPointOffset(int sampleNo,  DataTagged::getPointOffset(int sampleNo,
505                             int dataPointNo) const                             int dataPointNo) const
506  {  {
507    int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);    int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
508    DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));    DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));
509    DataArrayView::ValueType::size_type offset=m_defaultValueOffset;    DataTypes::ValueType::size_type offset=m_defaultValueOffset;
510    if (pos!=m_offsetLookup.end()) {    if (pos!=m_offsetLookup.end()) {
511      offset=pos->second;      offset=pos->second;
512    }    }
513    return offset;    return offset;
514  }  }
515    
516  DataArrayView  DataTypes::ValueType::size_type
517  DataTagged::getDataPointByTag(int tag) const  DataTagged::getPointOffset(int sampleNo,
518                               int dataPointNo)
519  {  {
520    DataMapType::const_iterator pos(m_offsetLookup.find(tag));    int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
521    DataArrayView::ValueType::size_type offset=m_defaultValueOffset;    DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));
522      DataTypes::ValueType::size_type offset=m_defaultValueOffset;
523    if (pos!=m_offsetLookup.end()) {    if (pos!=m_offsetLookup.end()) {
524      offset=pos->second;      offset=pos->second;
525    }    }
526    DataArrayView temp(getPointDataView());    return offset;
   temp.setOffset(offset);  
   return temp;  
527  }  }
528    
529  DataArrayView  DataTypes::ValueType::size_type
530  DataTagged::getDataPoint(int sampleNo,  DataTagged::getOffsetForTag(int tag) const
                          int dataPointNo)  
531  {  {
532    EsysAssert(validSampleNo(sampleNo),"(getDataPoint) Invalid sampleNo: " << sampleNo);    DataMapType::const_iterator pos(m_offsetLookup.find(tag));
533    int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);    DataTypes::ValueType::size_type offset=m_defaultValueOffset;
534    return getDataPointByTag(tagKey);    if (pos!=m_offsetLookup.end()) {
535        offset=pos->second;
536      }
537      return offset;
538  }  }
539    
540  int  DataTypes::ValueType::const_reference
541  DataTagged::archiveData(ofstream& archiveFile,  DataTagged::getDataByTagRO(int tag, DataTypes::ValueType::size_type i) const
                         const DataArrayView::ValueType::size_type noValues) const  
542  {  {
543    return(m_data.archiveData(archiveFile, noValues));    DataMapType::const_iterator pos(m_offsetLookup.find(tag));
544      DataTypes::ValueType::size_type offset=m_defaultValueOffset;
545      if (pos!=m_offsetLookup.end()) {
546        offset=pos->second;
547      }
548      return m_data[offset+i];
549  }  }
550    
551  int  DataTypes::ValueType::reference
552  DataTagged::extractData(ifstream& archiveFile,  DataTagged::getDataByTagRW(int tag, DataTypes::ValueType::size_type i)
                         const DataArrayView::ValueType::size_type noValues)  
553  {  {
554    return(m_data.extractData(archiveFile, noValues));    CHECK_FOR_EX_WRITE
555      DataMapType::const_iterator pos(m_offsetLookup.find(tag));
556      DataTypes::ValueType::size_type offset=m_defaultValueOffset;
557      if (pos!=m_offsetLookup.end()) {
558        offset=pos->second;
559      }
560      return m_data[offset+i];
561  }  }
562    
563  void  void
564  DataTagged::symmetric(DataAbstract* ev)  DataTagged::symmetric(DataAbstract* ev)
565  {  {
# Line 427  DataTagged::symmetric(DataAbstract* ev) Line 570  DataTagged::symmetric(DataAbstract* ev)
570    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
571    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
572    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
573      ValueType& evVec=temp_ev->getVectorRW();
574      const ShapeType& evShape=temp_ev->getShape();
575    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
576        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
577        DataArrayView thisView=getDataPointByTag(i->first);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
578        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
579        DataArrayView::symmetric(thisView,0,evView,0);        DataMaths::symmetric(m_data,getShape(),offset,evVec, evShape, evoffset);
580    }    }
581    DataArrayView::symmetric(getDefaultValue(),0,temp_ev->getDefaultValue(),0);    DataMaths::symmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
582  }  }
583    
584    
585  void  void
586  DataTagged::nonsymmetric(DataAbstract* ev)  DataTagged::nonsymmetric(DataAbstract* ev)
587  {  {
# Line 445  DataTagged::nonsymmetric(DataAbstract* e Line 592  DataTagged::nonsymmetric(DataAbstract* e
592    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
593    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
594    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
595      ValueType& evVec=temp_ev->getVectorRW();
596      const ShapeType& evShape=temp_ev->getShape();
597    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
598        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
599        DataArrayView thisView=getDataPointByTag(i->first);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
600        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
601        DataArrayView::nonsymmetric(thisView,0,evView,0);        DataMaths::nonsymmetric(m_data,getShape(),offset,evVec, evShape, evoffset);
602    }    }
603    DataArrayView::nonsymmetric(getDefaultValue(),0,temp_ev->getDefaultValue(),0);    DataMaths::nonsymmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
604  }  }
605    
606    
607  void  void
608  DataTagged::trace(DataAbstract* ev, int axis_offset)  DataTagged::trace(DataAbstract* ev, int axis_offset)
609  {  {
# Line 463  DataTagged::trace(DataAbstract* ev, int Line 614  DataTagged::trace(DataAbstract* ev, int
614    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
615    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
616    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
617      ValueType& evVec=temp_ev->getVectorRW();
618      const ShapeType& evShape=temp_ev->getShape();
619    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
620        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
621        DataArrayView thisView=getDataPointByTag(i->first);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
622        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
623        DataArrayView::trace(thisView,0,evView,0, axis_offset);        DataMaths::trace(m_data,getShape(),offset,evVec, evShape, evoffset, axis_offset);
624    }    }
625    DataArrayView::trace(getDefaultValue(),0,temp_ev->getDefaultValue(),0,axis_offset);    DataMaths::trace(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
626  }  }
627    
628  void  void
# Line 482  DataTagged::transpose(DataAbstract* ev, Line 635  DataTagged::transpose(DataAbstract* ev,
635    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
636    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
637    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
638      ValueType& evVec=temp_ev->getVectorRW();
639      const ShapeType& evShape=temp_ev->getShape();
640    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
641        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
642        DataArrayView thisView=getDataPointByTag(i->first);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
643        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
644        DataArrayView::transpose(thisView,0,evView,0, axis_offset);        DataMaths::transpose(m_data,getShape(),offset,evVec, evShape, evoffset, axis_offset);
645    }    }
646    DataArrayView::transpose(getDefaultValue(),0,temp_ev->getDefaultValue(),0,axis_offset);    DataMaths::transpose(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
647  }  }
648    
649  void  void
# Line 501  DataTagged::swapaxes(DataAbstract* ev, i Line 656  DataTagged::swapaxes(DataAbstract* ev, i
656    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
657    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
658    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
659      ValueType& evVec=temp_ev->getVectorRW();
660      const ShapeType& evShape=temp_ev->getShape();
661    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
662        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
663        DataArrayView thisView=getDataPointByTag(i->first);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
664        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
665        DataArrayView::swapaxes(thisView,0,evView,0,axis0,axis1);        DataMaths::swapaxes(m_data,getShape(),offset,evVec, evShape, evoffset,axis0,axis1);
666    }    }
667    DataArrayView::swapaxes(getDefaultValue(),0,temp_ev->getDefaultValue(),0,axis0,axis1);    DataMaths::swapaxes(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis0,axis1);
668  }  }
669    
670  void  void
# Line 520  DataTagged::eigenvalues(DataAbstract* ev Line 677  DataTagged::eigenvalues(DataAbstract* ev
677    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
678    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
679    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
680      ValueType& evVec=temp_ev->getVectorRW();
681      const ShapeType& evShape=temp_ev->getShape();
682    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
683        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
684        DataArrayView thisView=getDataPointByTag(i->first);  //       DataArrayView thisView=getDataPointByTag(i->first);
685        DataArrayView evView=temp_ev->getDataPointByTag(i->first);  //       DataArrayView evView=temp_ev->getDataPointByTag(i->first);
686        DataArrayView::eigenvalues(thisView,0,evView,0);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
687          DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
688          DataMaths::eigenvalues(m_data,getShape(),offset,evVec, evShape, evoffset);
689    }    }
690    DataArrayView::eigenvalues(getDefaultValue(),0,temp_ev->getDefaultValue(),0);    DataMaths::eigenvalues(m_data,getShape(),getDefaultOffset(),evVec, evShape, temp_ev->getDefaultOffset());
691  }  }
692  void  void
693  DataTagged::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)  DataTagged::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
# Line 542  DataTagged::eigenvalues_and_eigenvectors Line 703  DataTagged::eigenvalues_and_eigenvectors
703    const DataTagged::DataMapType& thisLookup=getTagLookup();    const DataTagged::DataMapType& thisLookup=getTagLookup();
704    DataTagged::DataMapType::const_iterator i;    DataTagged::DataMapType::const_iterator i;
705    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();    DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
706      ValueType& evVec=temp_ev->getVectorRW();
707      const ShapeType& evShape=temp_ev->getShape();
708      ValueType& VVec=temp_V->getVectorRW();
709      const ShapeType& VShape=temp_V->getShape();
710    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {    for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
711        temp_ev->addTaggedValue(i->first,temp_ev->getDefaultValue());        temp_ev->addTag(i->first);
712        temp_V->addTaggedValue(i->first,temp_V->getDefaultValue());        temp_V->addTag(i->first);
713        DataArrayView thisView=getDataPointByTag(i->first);  /*      DataArrayView thisView=getDataPointByTag(i->first);
714        DataArrayView evView=temp_ev->getDataPointByTag(i->first);        DataArrayView evView=temp_ev->getDataPointByTag(i->first);
715        DataArrayView VView=temp_V->getDataPointByTag(i->first);        DataArrayView VView=temp_V->getDataPointByTag(i->first);*/
716        DataArrayView::eigenvalues_and_eigenvectors(thisView,0,evView,0,VView,0,tol);        DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
717          DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
718          DataTypes::ValueType::size_type Voffset=temp_V->getOffsetForTag(i->first);
719    /*      DataArrayView::eigenvalues_and_eigenvectors(thisView,0,evView,0,VView,0,tol);*/
720          DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),offset,evVec, evShape, evoffset,VVec,VShape,Voffset,tol);
721    
722    }    }
723    DataArrayView::eigenvalues_and_eigenvectors(getDefaultValue(),0,    DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),getDefaultOffset(),evVec, evShape,
724                                                temp_ev->getDefaultValue(),0,                        temp_ev->getDefaultOffset(),VVec,VShape,
725                                                temp_V->getDefaultValue(),0,                        temp_V->getDefaultOffset(), tol);
                                               tol);  
726    
727    
728  }  }
729    
730  void  void
731    DataTagged::matrixInverse(DataAbstract* out) const
732    {
733      DataTagged* temp=dynamic_cast<DataTagged*>(out);
734      if (temp==0)
735      {
736        throw DataException("Error - DataTagged::matrixInverse: casting to DataTagged failed (propably a programming error).");
737      }
738      if (getRank()!=2)
739      {
740        throw DataException("Error - DataExpanded::matrixInverse: input must be rank 2.");
741      }
742      const DataTagged::DataMapType& thisLookup=getTagLookup();
743      DataTagged::DataMapType::const_iterator i;
744      DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
745      ValueType& outVec=temp->getVectorRW();
746      const ShapeType& outShape=temp->getShape();
747      LapackInverseHelper h(getShape()[0]);
748      int err=0;
749      for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
750          temp->addTag(i->first);
751          DataTypes::ValueType::size_type inoffset=getOffsetForTag(i->first);
752          DataTypes::ValueType::size_type outoffset=temp->getOffsetForTag(i->first);
753    
754          err=DataMaths::matrix_inverse(m_data, getShape(), inoffset, outVec, outShape, outoffset, 1, h);
755          if (!err) break;
756      }
757      if (!err)
758      {
759          DataMaths::matrix_inverse(m_data, getShape(), getDefaultOffset(), outVec, outShape, temp->getDefaultOffset(), 1, h);
760      }
761      if (err)
762      {
763         DataMaths::matrixInverseError(err);    // throws exceptions
764      }
765    }
766    
767    void
768  DataTagged::setToZero(){  DataTagged::setToZero(){
769      DataArrayView::ValueType::size_type n=m_data.size();      CHECK_FOR_EX_WRITE
770        DataTypes::ValueType::size_type n=m_data.size();
771      for (int i=0; i<n ;++i) m_data[i]=0.;      for (int i=0; i<n ;++i) m_data[i]=0.;
772  }  }
773    
774  void  void
775  DataTagged::dump(const std::string fileName) const  DataTagged::dump(const std::string fileName) const
776  {  {
    #ifdef PASO_MPI  
    throw DataException("Error - DataTagged:: dump is not implemented for MPI yet.");  
    #endif  
777     #ifdef USE_NETCDF     #ifdef USE_NETCDF
778     const int ldims=DataArrayView::maxRank+1;     const int ldims=DataTypes::maxRank+1;
779     const NcDim* ncdims[ldims];     const NcDim* ncdims[ldims];
780     NcVar *var, *tags_var;     NcVar *var, *tags_var;
781     int rank = getPointDataView().getRank();     int rank = getRank();
782     int type=  getFunctionSpace().getTypeCode();     int type=  getFunctionSpace().getTypeCode();
783     int ndims =0;     int ndims =0;
784     long dims[ldims];     long dims[ldims];
785     const double* d_ptr=&(m_data[0]);     const double* d_ptr=&(m_data[0]);
786     DataArrayView::ShapeType shape = getPointDataView().getShape();     DataTypes::ShapeType shape = getShape();
787       int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
788       int mpi_num=getFunctionSpace().getDomain()->getMPISize();
789    #ifdef PASO_MPI
790       MPI_Status status;
791    #endif
792    
793    #ifdef PASO_MPI
794       /* Serialize NetCDF I/O */
795       if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81803, MPI_COMM_WORLD, &status);
796    #endif
797    
798     // netCDF error handler     // netCDF error handler
799     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
800     // Create the file.     // Create the file.
801     NcFile dataFile(fileName.c_str(), NcFile::Replace);     char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
802       NcFile dataFile(newFileName, NcFile::Replace);
803     // check if writing was successful     // check if writing was successful
804     if (!dataFile.is_valid())     if (!dataFile.is_valid())
805          throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");          throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");
# Line 598  DataTagged::dump(const std::string fileN Line 813  DataTagged::dump(const std::string fileN
813     if ( rank >0 ) {     if ( rank >0 ) {
814         dims[0]=shape[0];         dims[0]=shape[0];
815         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
816              throw DataException("Error - DataTagged:: appending ncdimsion 0 to netCDF file failed.");              throw DataException("Error - DataTagged:: appending ncdimension 0 to netCDF file failed.");
817     }     }
818     if ( rank >1 ) {     if ( rank >1 ) {
819         dims[1]=shape[1];         dims[1]=shape[1];
820         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
821              throw DataException("Error - DataTagged:: appending ncdimsion 1 to netCDF file failed.");              throw DataException("Error - DataTagged:: appending ncdimension 1 to netCDF file failed.");
822     }     }
823     if ( rank >2 ) {     if ( rank >2 ) {
824         dims[2]=shape[2];         dims[2]=shape[2];
825         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
826              throw DataException("Error - DataTagged:: appending ncdimsion 2 to netCDF file failed.");              throw DataException("Error - DataTagged:: appending ncdimension 2 to netCDF file failed.");
827     }     }
828     if ( rank >3 ) {     if ( rank >3 ) {
829         dims[3]=shape[3];         dims[3]=shape[3];
830         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
831              throw DataException("Error - DataTagged:: appending ncdimsion 3 to netCDF file failed.");              throw DataException("Error - DataTagged:: appending ncdimension 3 to netCDF file failed.");
832     }     }
833     const DataTagged::DataMapType& thisLookup=getTagLookup();     const DataTagged::DataMapType& thisLookup=getTagLookup();
834     DataTagged::DataMapType::const_iterator i;     DataTagged::DataMapType::const_iterator i;
835     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
836     int ntags=1;     int ntags=1;
837     for (i=thisLookup.begin();i!=thisLookupEnd;i++) ntags++;     for (i=thisLookup.begin();i!=thisLookupEnd;i++) ntags++;
838     int* tags =(int*) malloc(ntags*sizeof(int));     int* tags =(int*) esysUtils::malloc(ntags*sizeof(int));
839     int c=1;     int c=1;
840     tags[0]=-1;     tags[0]=-1;
841     for (i=thisLookup.begin();i!=thisLookupEnd;i++) tags[c++]=i->first;     for (i=thisLookup.begin();i!=thisLookupEnd;i++) tags[c++]=i->first;
842     dims[rank]=ntags;     dims[rank]=ntags;
843     if (! (ncdims[rank] = dataFile.add_dim("num_tags", dims[rank])) )     if (! (ncdims[rank] = dataFile.add_dim("num_tags", dims[rank])) )
844     {     {
845         free(tags);         esysUtils::free(tags);
846             throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");             throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");
847     }     }
848     if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )     if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )
849     {     {
850      free(tags);      esysUtils::free(tags);
851          throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");          throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");
852     }     }
853     if (! (tags_var->put(tags,dims[rank])) )     if (! (tags_var->put(tags,dims[rank])) )
854     {     {
855      free(tags);      esysUtils::free(tags);
856          throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");          throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");
857     }     }
858     if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )     if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
859     {     {
860      free(tags);      esysUtils::free(tags);
861          throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");          throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");
862     }     }
863     if (! (var->put(d_ptr,dims)) )     if (! (var->put(d_ptr,dims)) )
864     {     {
865      free(tags);      esysUtils::free(tags);
866          throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");          throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");
867     }     }
868    #ifdef PASO_MPI
869       if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81803, MPI_COMM_WORLD);
870    #endif
871     #else     #else
872     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.");
873     #endif     #endif
874  }  }
875    
876    DataTypes::ValueType&
877    DataTagged::getVectorRW()
878    {
879        CHECK_FOR_EX_WRITE
880        return m_data;
881    }
882    
883    const DataTypes::ValueType&
884    DataTagged::getVectorRO() const
885    {
886        return m_data;
887    }
888    
889  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1455  
changed lines
  Added in v.2766

  ViewVC Help
Powered by ViewVC 1.1.26