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

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

Legend:
Removed from v.119  
changed lines
  Added in v.1834

  ViewVC Help
Powered by ViewVC 1.1.26