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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2005 - (hide annotations)
Mon Nov 10 01:21:39 2008 UTC (11 years, 3 months ago) by jfenwick
File size: 35373 byte(s)
Bringing all changes across from schroedinger.
(Note this does not mean development is done, just that it will happen
on the trunk for now).
If anyone notices any problems please contact me.


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 ksteube 1827 #include "Data.h"
16 jgs 474 #include "DataTagged.h"
17 phornby 1628 #include "esysUtils/esys_malloc.h"
18 jgs 496
19 jgs 474 #include "DataConstant.h"
20     #include "DataException.h"
21 gross 1023 #ifdef USE_NETCDF
22 ksteube 1312 #include <netcdfcpp.h>
23 gross 1023 #endif
24 ksteube 1827 #ifdef PASO_MPI
25     #include <mpi.h>
26     #endif
27 jfenwick 1796 #include "DataMaths.h"
28 jgs 82
29     using namespace std;
30    
31     namespace escript {
32    
33 jgs 121 DataTagged::DataTagged()
34 jfenwick 2005 : parent(FunctionSpace(),DataTypes::scalarShape)
35 jgs 82 {
36 jgs 509 // default constructor
37    
38 jgs 82 // create a scalar default value
39 jgs 151 m_data.resize(1,0.,1);
40 jfenwick 1796 /* DataArrayView temp(m_data,DataTypes::ShapeType());
41     setPointDataView(temp);*/
42 jgs 82 }
43    
44 jfenwick 1796 // 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 jgs 509
68 jgs 119 DataTagged::DataTagged(const FunctionSpace& what,
69 jfenwick 1796 const DataTypes::ShapeType &shape,
70 jgs 119 const int tags[],
71 jgs 500 const ValueType& data)
72 jfenwick 2005 : parent(what,shape)
73 jgs 119 {
74 jgs 509 // alternative constructor
75     // not unit_tested tested yet
76 jfenwick 1796 // It is not explicitly unit tested yet, but it is called from DataFactory
77 jgs 509
78 jfenwick 1802 if (!what.canTag())
79     {
80     throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
81     }
82 jgs 121 // copy the data
83 jgs 119 m_data=data;
84 jgs 121
85 jgs 119 // create the view of the data
86 jfenwick 1796 // DataArrayView tempView(m_data,shape);
87     // setPointDataView(tempView);
88 jgs 121
89 jfenwick 1796 // 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 jgs 119 // create the tag lookup map
95 jfenwick 1796 // 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 jgs 119 }
100     }
101    
102 woo409 757 DataTagged::DataTagged(const FunctionSpace& what,
103 jfenwick 1796 const DataTypes::ShapeType &shape,
104 woo409 757 const TagListType& tags,
105     const ValueType& data)
106 jfenwick 2005 : parent(what,shape)
107 woo409 757 {
108     // alternative constructor
109    
110 jfenwick 1802 if (!what.canTag())
111     {
112     throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
113     }
114    
115 woo409 757 // copy the data
116     m_data=data;
117    
118     // create the view of the data
119 jfenwick 1796 // DataArrayView tempView(m_data,shape);
120     // setPointDataView(tempView);
121 woo409 757
122     // create the tag lookup map
123 jfenwick 1796
124     // 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 woo409 757 }
137 jfenwick 1796
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 woo409 757 }
145    
146    
147 jgs 102 DataTagged::DataTagged(const DataTagged& other)
148 jfenwick 2005 : parent(other.getFunctionSpace(),other.getShape()),
149 jfenwick 1946 m_offsetLookup(other.m_offsetLookup),
150     m_data(other.m_data)
151 jgs 82 {
152 jgs 509 // copy constructor
153    
154 jgs 119 // create the data view
155 jfenwick 1796 // DataArrayView temp(m_data,other.getPointDataView().getShape());
156     // setPointDataView(temp);
157 jgs 82 }
158    
159 jgs 102 DataTagged::DataTagged(const DataConstant& other)
160 jfenwick 2005 : parent(other.getFunctionSpace(),other.getShape())
161 jgs 82 {
162 jgs 509 // copy constructor
163    
164 jfenwick 1802 if (!other.getFunctionSpace().canTag())
165     {
166     throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
167     }
168    
169 jgs 121 // fill the default value with the constant value item from "other"
170 jfenwick 1796 // const DataArrayView& value=other.getPointDataView();
171     int len = other.getNoValues();
172 jgs 151 m_data.resize(len,0.,len);
173 jfenwick 1796 for (int i=0; i<len; i++) {
174     m_data[i]=other.getVector()[i];
175 jgs 121 }
176    
177 jfenwick 1796 // // create the data view
178     // DataArrayView temp(m_data,value.getShape());
179     // setPointDataView(temp);
180 jgs 82 }
181    
182 jfenwick 1796
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 jfenwick 2005 : parent(what,shape)
189 jfenwick 1796 {
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 jfenwick 1802 if (!what.canTag())
198     {
199     throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
200     }
201    
202 jfenwick 1796 if (tagsource!=0)
203     {
204 jfenwick 1946 m_data.resize(defaultvalue.size(),0.); // since this is tagged data, we should have blocksize=1
205 jfenwick 1796
206     DataTagged::DataMapType::const_iterator i;
207     for (i=tagsource->getTagLookup().begin();i!=tagsource->getTagLookup().end();i++) {
208     addTag(i->first);
209     }
210     }
211     else
212     {
213     m_data.resize(defaultvalue.size());
214     }
215    
216     // need to set the default value ....
217     for (int i=0; i<defaultvalue.size(); i++) {
218     m_data[i]=defaultvalue[i];
219     }
220     }
221    
222 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 jfenwick 2005 : parent(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 phornby 1834 for (size_t i=0;i<values.size();++i)
443 jfenwick 1796 {
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 1981 unsigned int n=getNoValues();
457     unsigned 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 2005 DataTypes::ValueType::size_type
632     DataTagged::getPointOffset(int sampleNo,
633     int dataPointNo)
634     {
635     int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
636     DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));
637     DataTypes::ValueType::size_type offset=m_defaultValueOffset;
638     if (pos!=m_offsetLookup.end()) {
639     offset=pos->second;
640     }
641     return offset;
642     }
643    
644 jfenwick 1796 // DataArrayView
645     // DataTagged::getDataPointByTag(int tag) const
646     // {
647     // DataMapType::const_iterator pos(m_offsetLookup.find(tag));
648     // DataTypes::ValueType::size_type offset=m_defaultValueOffset;
649     // if (pos!=m_offsetLookup.end()) {
650     // offset=pos->second;
651     // }
652     // DataArrayView temp(getPointDataView());
653     // temp.setOffset(offset);
654     // return temp;
655     // }
656     //
657    
658    
659     DataTypes::ValueType::size_type
660     DataTagged::getOffsetForTag(int tag) const
661 jgs 82 {
662 jgs 496 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
663 jfenwick 1796 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
664 jgs 82 if (pos!=m_offsetLookup.end()) {
665     offset=pos->second;
666     }
667 jfenwick 1796 return offset;
668 jgs 82 }
669    
670 jfenwick 1796 DataTypes::ValueType::const_reference
671     DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i) const
672 jgs 82 {
673 jfenwick 1796 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
674     DataTypes::ValueType::size_type offset=m_defaultValueOffset;
675     if (pos!=m_offsetLookup.end()) {
676     offset=pos->second;
677     }
678     return m_data[offset+i];
679     /* DataArrayView temp(getPointDataView());
680     temp.setOffset(offset);
681     return temp.getData()[offset+i];*/
682 jgs 82 }
683    
684 jfenwick 1796
685     DataTypes::ValueType::reference
686     DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i)
687 jgs 123 {
688 jfenwick 1796 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
689     DataTypes::ValueType::size_type offset=m_defaultValueOffset;
690     if (pos!=m_offsetLookup.end()) {
691     offset=pos->second;
692     }
693     return m_data[offset+i];
694     /* DataArrayView temp(getPointDataView());
695     temp.setOffset(offset);
696     return temp.getData()[offset+i];*/
697 jgs 123 }
698    
699 jfenwick 1796
700    
701    
702    
703    
704     // DataArrayView
705     // DataTagged::getDataPoint(int sampleNo,
706     // int dataPointNo)
707     // {
708     // EsysAssert(validSampleNo(sampleNo),"(getDataPoint) Invalid sampleNo: " << sampleNo);
709     // int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
710     // return getDataPointByTag(tagKey);
711     // }
712    
713    
714 gross 580 void
715 ksteube 775 DataTagged::symmetric(DataAbstract* ev)
716     {
717     DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
718     if (temp_ev==0) {
719     throw DataException("Error - DataTagged::symmetric casting to DataTagged failed (probably a programming error).");
720     }
721     const DataTagged::DataMapType& thisLookup=getTagLookup();
722     DataTagged::DataMapType::const_iterator i;
723     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
724 jfenwick 1796 ValueType& evVec=temp_ev->getVector();
725     const ShapeType& evShape=temp_ev->getShape();
726 ksteube 775 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
727 jfenwick 1796 temp_ev->addTag(i->first);
728     DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
729     // DataArrayView thisView=getDataPointByTag(i->first);
730     // DataArrayView evView=temp_ev->getDataPointByTag(i->first);
731     DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
732    
733     // DataArrayView::symmetric(thisView,0,evView,0);
734     DataMaths::symmetric(m_data,getShape(),offset,evVec, evShape, evoffset);
735 ksteube 775 }
736 jfenwick 1796 // symmetric(m_data,getShape(),getDefaultOffset(),
737     DataMaths::symmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
738 ksteube 775 }
739 jfenwick 1796
740    
741 ksteube 775 void
742     DataTagged::nonsymmetric(DataAbstract* ev)
743     {
744     DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
745     if (temp_ev==0) {
746     throw DataException("Error - DataTagged::nonsymmetric casting to DataTagged failed (probably a programming error).");
747     }
748     const DataTagged::DataMapType& thisLookup=getTagLookup();
749     DataTagged::DataMapType::const_iterator i;
750     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
751 jfenwick 1796 ValueType& evVec=temp_ev->getVector();
752     const ShapeType& evShape=temp_ev->getShape();
753 ksteube 775 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
754 jfenwick 1796 temp_ev->addTag(i->first);
755     /* DataArrayView thisView=getDataPointByTag(i->first);
756     DataArrayView evView=temp_ev->getDataPointByTag(i->first);*/
757     DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
758     DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
759     DataMaths::nonsymmetric(m_data,getShape(),offset,evVec, evShape, evoffset);
760 ksteube 775 }
761 jfenwick 1796 DataMaths::nonsymmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
762 ksteube 775 }
763 jfenwick 1796
764    
765 ksteube 775 void
766 gross 800 DataTagged::trace(DataAbstract* ev, int axis_offset)
767 ksteube 775 {
768     DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
769     if (temp_ev==0) {
770 gross 800 throw DataException("Error - DataTagged::trace casting to DataTagged failed (probably a programming error).");
771 ksteube 775 }
772     const DataTagged::DataMapType& thisLookup=getTagLookup();
773     DataTagged::DataMapType::const_iterator i;
774     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
775 jfenwick 1796 ValueType& evVec=temp_ev->getVector();
776     const ShapeType& evShape=temp_ev->getShape();
777 ksteube 775 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
778 jfenwick 1796 temp_ev->addTag(i->first);
779     // DataArrayView thisView=getDataPointByTag(i->first);
780     // DataArrayView evView=temp_ev->getDataPointByTag(i->first);
781     DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
782     DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
783     DataMaths::trace(m_data,getShape(),offset,evVec, evShape, evoffset, axis_offset);
784 ksteube 775 }
785 jfenwick 1796 DataMaths::trace(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
786 ksteube 775 }
787 gross 800
788 ksteube 775 void
789     DataTagged::transpose(DataAbstract* ev, int axis_offset)
790     {
791     DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
792     if (temp_ev==0) {
793     throw DataException("Error - DataTagged::transpose casting to DataTagged failed (probably a programming error).");
794     }
795     const DataTagged::DataMapType& thisLookup=getTagLookup();
796     DataTagged::DataMapType::const_iterator i;
797     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
798 jfenwick 1796 ValueType& evVec=temp_ev->getVector();
799     const ShapeType& evShape=temp_ev->getShape();
800 ksteube 775 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
801 jfenwick 1796 temp_ev->addTag(i->first);
802     // DataArrayView thisView=getDataPointByTag(i->first);
803     // DataArrayView evView=temp_ev->getDataPointByTag(i->first);
804     DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
805     DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
806     DataMaths::transpose(m_data,getShape(),offset,evVec, evShape, evoffset, axis_offset);
807 ksteube 775 }
808 jfenwick 1796 DataMaths::transpose(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
809 ksteube 775 }
810 gross 800
811 ksteube 775 void
812 gross 804 DataTagged::swapaxes(DataAbstract* ev, int axis0, int axis1)
813 gross 800 {
814     DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
815     if (temp_ev==0) {
816 gross 804 throw DataException("Error - DataTagged::swapaxes casting to DataTagged failed (probably a programming error).");
817 gross 800 }
818     const DataTagged::DataMapType& thisLookup=getTagLookup();
819     DataTagged::DataMapType::const_iterator i;
820     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
821 jfenwick 1796 ValueType& evVec=temp_ev->getVector();
822     const ShapeType& evShape=temp_ev->getShape();
823 gross 800 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
824 jfenwick 1796 temp_ev->addTag(i->first);
825     /* DataArrayView thisView=getDataPointByTag(i->first);
826     DataArrayView evView=temp_ev->getDataPointByTag(i->first);*/
827     DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
828     DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
829     DataMaths::swapaxes(m_data,getShape(),offset,evVec, evShape, evoffset,axis0,axis1);
830 gross 800 }
831 jfenwick 1796 DataMaths::swapaxes(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis0,axis1);
832 gross 800 }
833    
834     void
835 gross 580 DataTagged::eigenvalues(DataAbstract* ev)
836     {
837     DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
838     if (temp_ev==0) {
839     throw DataException("Error - DataTagged::eigenvalues casting to DataTagged failed (propably a programming error).");
840     }
841     const DataTagged::DataMapType& thisLookup=getTagLookup();
842     DataTagged::DataMapType::const_iterator i;
843     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
844 jfenwick 1796 ValueType& evVec=temp_ev->getVector();
845     const ShapeType& evShape=temp_ev->getShape();
846 gross 580 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
847 jfenwick 1796 temp_ev->addTag(i->first);
848     // DataArrayView thisView=getDataPointByTag(i->first);
849     // DataArrayView evView=temp_ev->getDataPointByTag(i->first);
850     DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
851     DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
852     DataMaths::eigenvalues(m_data,getShape(),offset,evVec, evShape, evoffset);
853 gross 580 }
854 jfenwick 1796 DataMaths::eigenvalues(m_data,getShape(),getDefaultOffset(),evVec, evShape, temp_ev->getDefaultOffset());
855 gross 580 }
856     void
857     DataTagged::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
858     {
859     DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
860     if (temp_ev==0) {
861     throw DataException("Error - DataTagged::eigenvalues_and_eigenvectors casting to DataTagged failed (propably a programming error).");
862     }
863     DataTagged* temp_V=dynamic_cast<DataTagged*>(V);
864     if (temp_V==0) {
865     throw DataException("Error - DataTagged::eigenvalues_and_eigenvectors casting to DataTagged failed (propably a programming error).");
866     }
867     const DataTagged::DataMapType& thisLookup=getTagLookup();
868     DataTagged::DataMapType::const_iterator i;
869     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
870 jfenwick 1796 ValueType& evVec=temp_ev->getVector();
871     const ShapeType& evShape=temp_ev->getShape();
872     ValueType& VVec=temp_V->getVector();
873     const ShapeType& VShape=temp_V->getShape();
874 gross 580 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
875 jfenwick 1796 temp_ev->addTag(i->first);
876     temp_V->addTag(i->first);
877     /* DataArrayView thisView=getDataPointByTag(i->first);
878 gross 580 DataArrayView evView=temp_ev->getDataPointByTag(i->first);
879 jfenwick 1796 DataArrayView VView=temp_V->getDataPointByTag(i->first);*/
880     DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
881     DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
882     DataTypes::ValueType::size_type Voffset=temp_V->getOffsetForTag(i->first);
883     /* DataArrayView::eigenvalues_and_eigenvectors(thisView,0,evView,0,VView,0,tol);*/
884     DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),offset,evVec, evShape, evoffset,VVec,VShape,Voffset,tol);
885    
886 gross 580 }
887 jfenwick 1796 DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),getDefaultOffset(),evVec, evShape,
888     temp_ev->getDefaultOffset(),VVec,VShape,
889     temp_V->getDefaultOffset(), tol);
890 jgs 123
891 gross 580
892     }
893    
894 gross 950 void
895 gross 1118 DataTagged::setToZero(){
896 jfenwick 1796 DataTypes::ValueType::size_type n=m_data.size();
897 gross 1118 for (int i=0; i<n ;++i) m_data[i]=0.;
898     }
899    
900     void
901 gross 950 DataTagged::dump(const std::string fileName) const
902     {
903 gross 1023 #ifdef USE_NETCDF
904 jfenwick 1796 const int ldims=DataTypes::maxRank+1;
905 gross 983 const NcDim* ncdims[ldims];
906     NcVar *var, *tags_var;
907 jfenwick 1796 int rank = getRank();
908 gross 983 int type= getFunctionSpace().getTypeCode();
909     int ndims =0;
910     long dims[ldims];
911 gross 1141 const double* d_ptr=&(m_data[0]);
912 jfenwick 1796 DataTypes::ShapeType shape = getShape();
913 jfenwick 1872 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
914     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
915 ksteube 1827 #ifdef PASO_MPI
916     MPI_Status status;
917     #endif
918 gross 983
919 ksteube 1827 #ifdef PASO_MPI
920     /* Serialize NetCDF I/O */
921     if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81803, MPI_COMM_WORLD, &status);
922     #endif
923    
924 gross 983 // netCDF error handler
925     NcError err(NcError::verbose_nonfatal);
926     // Create the file.
927 ksteube 1827 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
928     NcFile dataFile(newFileName, NcFile::Replace);
929 gross 983 // check if writing was successful
930     if (!dataFile.is_valid())
931     throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");
932 gross 1141 if (!dataFile.add_att("type_id",1) )
933 gross 983 throw DataException("Error - DataTagged:: appending data type to netCDF file failed.");
934     if (!dataFile.add_att("rank",rank) )
935     throw DataException("Error - DataTagged:: appending rank attribute to netCDF file failed.");
936     if (!dataFile.add_att("function_space_type",type))
937     throw DataException("Error - DataTagged:: appending function space attribute to netCDF file failed.");
938     ndims=rank+1;
939     if ( rank >0 ) {
940     dims[0]=shape[0];
941     if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
942 ksteube 1800 throw DataException("Error - DataTagged:: appending ncdimension 0 to netCDF file failed.");
943 gross 983 }
944     if ( rank >1 ) {
945     dims[1]=shape[1];
946     if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
947 ksteube 1800 throw DataException("Error - DataTagged:: appending ncdimension 1 to netCDF file failed.");
948 gross 983 }
949     if ( rank >2 ) {
950     dims[2]=shape[2];
951     if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
952 ksteube 1800 throw DataException("Error - DataTagged:: appending ncdimension 2 to netCDF file failed.");
953 gross 983 }
954     if ( rank >3 ) {
955     dims[3]=shape[3];
956     if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
957 ksteube 1800 throw DataException("Error - DataTagged:: appending ncdimension 3 to netCDF file failed.");
958 gross 983 }
959     const DataTagged::DataMapType& thisLookup=getTagLookup();
960     DataTagged::DataMapType::const_iterator i;
961     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
962     int ntags=1;
963     for (i=thisLookup.begin();i!=thisLookupEnd;i++) ntags++;
964 phornby 1628 int* tags =(int*) esysUtils::malloc(ntags*sizeof(int));
965 phornby 1020 int c=1;
966 gross 983 tags[0]=-1;
967     for (i=thisLookup.begin();i!=thisLookupEnd;i++) tags[c++]=i->first;
968     dims[rank]=ntags;
969     if (! (ncdims[rank] = dataFile.add_dim("num_tags", dims[rank])) )
970 phornby 1020 {
971 phornby 1628 esysUtils::free(tags);
972 phornby 1020 throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");
973     }
974 gross 983 if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )
975 phornby 1020 {
976 phornby 1628 esysUtils::free(tags);
977 gross 983 throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");
978 phornby 1020 }
979 gross 983 if (! (tags_var->put(tags,dims[rank])) )
980 phornby 1020 {
981 phornby 1628 esysUtils::free(tags);
982 gross 983 throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");
983 phornby 1020 }
984 gross 983 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
985 phornby 1020 {
986 phornby 1628 esysUtils::free(tags);
987 gross 983 throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");
988 phornby 1020 }
989 gross 1141 if (! (var->put(d_ptr,dims)) )
990 phornby 1020 {
991 phornby 1628 esysUtils::free(tags);
992 gross 983 throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");
993 phornby 1020 }
994 ksteube 1827 #ifdef PASO_MPI
995     if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81803, MPI_COMM_WORLD);
996     #endif
997 gross 1023 #else
998     throw DataException("Error - DataTagged:: dump is not configured with netCDF. Please contact your installation manager.");
999     #endif
1000 gross 950 }
1001 jfenwick 1796
1002     DataTypes::ValueType&
1003     DataTagged::getVector()
1004     {
1005     return m_data;
1006     }
1007    
1008     const DataTypes::ValueType&
1009     DataTagged::getVector() const
1010     {
1011     return m_data;
1012     }
1013    
1014 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