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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 1 month ago) by ksteube
File size: 34650 byte(s)
Copyright updated in all files

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26