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

Legend:
Removed from v.82  
changed lines
  Added in v.2742

  ViewVC Help
Powered by ViewVC 1.1.26