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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1800 - (hide annotations)
Thu Sep 18 05:28:20 2008 UTC (11 years, 1 month ago) by ksteube
File size: 34138 byte(s)
Serialized parallel I/O when writing mesh or data to NetCDF file on multiple MPI processors.
Added domain method getMPIComm() to complement getMPISize() and getMPIRank().

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