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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1981 - (hide annotations)
Thu Nov 6 05:27:33 2008 UTC (11 years, 1 month ago) by jfenwick
File size: 35023 byte(s)
More warning removal.

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 1796 : DataAbstract(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 1796 : DataAbstract(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 1796 : DataAbstract(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 1796 : DataAbstract(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 1796 : DataAbstract(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     : DataAbstract(what,shape)
189     {
190     // This constructor has not been unit tested yet
191    
192     if (defaultvalue.size()!=DataTypes::noValues(shape)) {
193     throw DataException("Programming error - defaultvalue does not match supplied shape.");
194     }
195    
196    
197 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     : 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 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 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 USE_NETCDF
891 jfenwick 1796 const int ldims=DataTypes::maxRank+1;
892 gross 983 const NcDim* ncdims[ldims];
893     NcVar *var, *tags_var;
894 jfenwick 1796 int rank = getRank();
895 gross 983 int type= getFunctionSpace().getTypeCode();
896     int ndims =0;
897     long dims[ldims];
898 gross 1141 const double* d_ptr=&(m_data[0]);
899 jfenwick 1796 DataTypes::ShapeType shape = getShape();
900 jfenwick 1872 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
901     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
902 ksteube 1827 #ifdef PASO_MPI
903     MPI_Status status;
904     #endif
905 gross 983
906 ksteube 1827 #ifdef PASO_MPI
907     /* Serialize NetCDF I/O */
908     if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81803, MPI_COMM_WORLD, &status);
909     #endif
910    
911 gross 983 // netCDF error handler
912     NcError err(NcError::verbose_nonfatal);
913     // Create the file.
914 ksteube 1827 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
915     NcFile dataFile(newFileName, NcFile::Replace);
916 gross 983 // check if writing was successful
917     if (!dataFile.is_valid())
918     throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");
919 gross 1141 if (!dataFile.add_att("type_id",1) )
920 gross 983 throw DataException("Error - DataTagged:: appending data type to netCDF file failed.");
921     if (!dataFile.add_att("rank",rank) )
922     throw DataException("Error - DataTagged:: appending rank attribute to netCDF file failed.");
923     if (!dataFile.add_att("function_space_type",type))
924     throw DataException("Error - DataTagged:: appending function space attribute to netCDF file failed.");
925     ndims=rank+1;
926     if ( rank >0 ) {
927     dims[0]=shape[0];
928     if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
929 ksteube 1800 throw DataException("Error - DataTagged:: appending ncdimension 0 to netCDF file failed.");
930 gross 983 }
931     if ( rank >1 ) {
932     dims[1]=shape[1];
933     if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
934 ksteube 1800 throw DataException("Error - DataTagged:: appending ncdimension 1 to netCDF file failed.");
935 gross 983 }
936     if ( rank >2 ) {
937     dims[2]=shape[2];
938     if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
939 ksteube 1800 throw DataException("Error - DataTagged:: appending ncdimension 2 to netCDF file failed.");
940 gross 983 }
941     if ( rank >3 ) {
942     dims[3]=shape[3];
943     if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
944 ksteube 1800 throw DataException("Error - DataTagged:: appending ncdimension 3 to netCDF file failed.");
945 gross 983 }
946     const DataTagged::DataMapType& thisLookup=getTagLookup();
947     DataTagged::DataMapType::const_iterator i;
948     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
949     int ntags=1;
950     for (i=thisLookup.begin();i!=thisLookupEnd;i++) ntags++;
951 phornby 1628 int* tags =(int*) esysUtils::malloc(ntags*sizeof(int));
952 phornby 1020 int c=1;
953 gross 983 tags[0]=-1;
954     for (i=thisLookup.begin();i!=thisLookupEnd;i++) tags[c++]=i->first;
955     dims[rank]=ntags;
956     if (! (ncdims[rank] = dataFile.add_dim("num_tags", dims[rank])) )
957 phornby 1020 {
958 phornby 1628 esysUtils::free(tags);
959 phornby 1020 throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");
960     }
961 gross 983 if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )
962 phornby 1020 {
963 phornby 1628 esysUtils::free(tags);
964 gross 983 throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");
965 phornby 1020 }
966 gross 983 if (! (tags_var->put(tags,dims[rank])) )
967 phornby 1020 {
968 phornby 1628 esysUtils::free(tags);
969 gross 983 throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");
970 phornby 1020 }
971 gross 983 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
972 phornby 1020 {
973 phornby 1628 esysUtils::free(tags);
974 gross 983 throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");
975 phornby 1020 }
976 gross 1141 if (! (var->put(d_ptr,dims)) )
977 phornby 1020 {
978 phornby 1628 esysUtils::free(tags);
979 gross 983 throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");
980 phornby 1020 }
981 ksteube 1827 #ifdef PASO_MPI
982     if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81803, MPI_COMM_WORLD);
983     #endif
984 gross 1023 #else
985     throw DataException("Error - DataTagged:: dump is not configured with netCDF. Please contact your installation manager.");
986     #endif
987 gross 950 }
988 jfenwick 1796
989     DataTypes::ValueType&
990     DataTagged::getVector()
991     {
992     return m_data;
993     }
994    
995     const DataTypes::ValueType&
996     DataTagged::getVector() const
997     {
998     return m_data;
999     }
1000    
1001 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