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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1872 - (hide annotations)
Mon Oct 13 00:18:55 2008 UTC (11 years, 2 months ago) by jfenwick
File size: 35110 byte(s)
Closing the moreshared branch

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