1 |
|
2 |
/******************************************************* |
3 |
* |
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 |
|
14 |
|
15 |
#include "Data.h" |
16 |
#include "DataTagged.h" |
17 |
#include "esysUtils/esys_malloc.h" |
18 |
|
19 |
#include "DataConstant.h" |
20 |
#include "DataException.h" |
21 |
#ifdef USE_NETCDF |
22 |
#include <netcdfcpp.h> |
23 |
#endif |
24 |
#ifdef PASO_MPI |
25 |
#include <mpi.h> |
26 |
#endif |
27 |
#include "DataMaths.h" |
28 |
|
29 |
using namespace std; |
30 |
|
31 |
namespace escript { |
32 |
|
33 |
DataTagged::DataTagged() |
34 |
: DataAbstract(FunctionSpace(),DataTypes::scalarShape) |
35 |
{ |
36 |
// default constructor |
37 |
|
38 |
// create a scalar default value |
39 |
m_data.resize(1,0.,1); |
40 |
/* DataArrayView temp(m_data,DataTypes::ShapeType()); |
41 |
setPointDataView(temp);*/ |
42 |
} |
43 |
|
44 |
// 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 |
|
68 |
DataTagged::DataTagged(const FunctionSpace& what, |
69 |
const DataTypes::ShapeType &shape, |
70 |
const int tags[], |
71 |
const ValueType& data) |
72 |
: DataAbstract(what,shape) |
73 |
{ |
74 |
// alternative constructor |
75 |
// not unit_tested tested yet |
76 |
// It is not explicitly unit tested yet, but it is called from DataFactory |
77 |
|
78 |
if (!what.canTag()) |
79 |
{ |
80 |
throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace."); |
81 |
} |
82 |
// copy the data |
83 |
m_data=data; |
84 |
|
85 |
// create the view of the data |
86 |
// DataArrayView tempView(m_data,shape); |
87 |
// setPointDataView(tempView); |
88 |
|
89 |
// 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 |
// create the tag lookup map |
95 |
// 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 |
} |
100 |
} |
101 |
|
102 |
DataTagged::DataTagged(const FunctionSpace& what, |
103 |
const DataTypes::ShapeType &shape, |
104 |
const TagListType& tags, |
105 |
const ValueType& data) |
106 |
: DataAbstract(what,shape) |
107 |
{ |
108 |
// alternative constructor |
109 |
|
110 |
if (!what.canTag()) |
111 |
{ |
112 |
throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace."); |
113 |
} |
114 |
|
115 |
// copy the data |
116 |
m_data=data; |
117 |
|
118 |
// create the view of the data |
119 |
// DataArrayView tempView(m_data,shape); |
120 |
// setPointDataView(tempView); |
121 |
|
122 |
// create the tag lookup map |
123 |
|
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 |
} |
137 |
|
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 |
} |
145 |
|
146 |
|
147 |
DataTagged::DataTagged(const DataTagged& other) |
148 |
: DataAbstract(other.getFunctionSpace(),other.getShape()), |
149 |
m_offsetLookup(other.m_offsetLookup), |
150 |
m_data(other.m_data) |
151 |
{ |
152 |
// copy constructor |
153 |
|
154 |
// create the data view |
155 |
// DataArrayView temp(m_data,other.getPointDataView().getShape()); |
156 |
// setPointDataView(temp); |
157 |
} |
158 |
|
159 |
DataTagged::DataTagged(const DataConstant& other) |
160 |
: DataAbstract(other.getFunctionSpace(),other.getShape()) |
161 |
{ |
162 |
// copy constructor |
163 |
|
164 |
if (!other.getFunctionSpace().canTag()) |
165 |
{ |
166 |
throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace."); |
167 |
} |
168 |
|
169 |
// fill the default value with the constant value item from "other" |
170 |
// const DataArrayView& value=other.getPointDataView(); |
171 |
int len = other.getNoValues(); |
172 |
m_data.resize(len,0.,len); |
173 |
for (int i=0; i<len; i++) { |
174 |
m_data[i]=other.getVector()[i]; |
175 |
} |
176 |
|
177 |
// // create the data view |
178 |
// DataArrayView temp(m_data,value.getShape()); |
179 |
// setPointDataView(temp); |
180 |
} |
181 |
|
182 |
|
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 |
if (!what.canTag()) |
198 |
{ |
199 |
throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace."); |
200 |
} |
201 |
|
202 |
if (tagsource!=0) |
203 |
{ |
204 |
m_data.resize(defaultvalue.size(),0.); // since this is tagged data, we should have blocksize=1 |
205 |
|
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 |
DataAbstract* |
223 |
DataTagged::deepCopy() |
224 |
{ |
225 |
return new DataTagged(*this); |
226 |
} |
227 |
|
228 |
DataAbstract* |
229 |
DataTagged::getSlice(const DataTypes::RegionType& region) const |
230 |
{ |
231 |
return new DataTagged(*this, region); |
232 |
} |
233 |
|
234 |
DataTagged::DataTagged(const DataTagged& other, |
235 |
const DataTypes::RegionType& region) |
236 |
: DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region)) |
237 |
{ |
238 |
// slice constructor |
239 |
|
240 |
// get the shape of the slice to copy from other |
241 |
DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region)); |
242 |
DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region); |
243 |
|
244 |
// allocate enough space in this for all values |
245 |
// (need to add one to allow for the default value) |
246 |
int len = DataTypes::noValues(regionShape)*(other.m_offsetLookup.size()+1); |
247 |
m_data.resize(len,0.0,len); |
248 |
|
249 |
// create the data view |
250 |
// DataArrayView temp(m_data,regionShape); |
251 |
// setPointDataView(temp); |
252 |
|
253 |
// copy the default value from other to this |
254 |
// 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 |
|
259 |
// loop through the tag values copying these |
260 |
DataMapType::const_iterator pos; |
261 |
DataTypes::ValueType::size_type tagOffset=getNoValues(); |
262 |
for (pos=other.m_offsetLookup.begin();pos!=other.m_offsetLookup.end();pos++){ |
263 |
// getPointDataView().copySlice(tagOffset,other.getPointDataView(),pos->second,regionLoopRange); |
264 |
DataTypes::copySlice(m_data,getShape(),tagOffset,otherData, otherShape, pos->second, regionLoopRange); |
265 |
m_offsetLookup.insert(DataMapType::value_type(pos->first,tagOffset)); |
266 |
tagOffset+=getNoValues(); |
267 |
} |
268 |
} |
269 |
|
270 |
void |
271 |
DataTagged::setSlice(const DataAbstract* other, |
272 |
const DataTypes::RegionType& region) |
273 |
{ |
274 |
|
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 |
throw DataException("Programming error - casting to DataTagged."); |
280 |
} |
281 |
|
282 |
// determine shape of the specified region |
283 |
DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region)); |
284 |
|
285 |
// modify region specification as needed to match rank of this object |
286 |
DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region); |
287 |
|
288 |
// ensure rank/shape of this object is compatible with specified region |
289 |
if (getRank()!=region.size()) { |
290 |
throw DataException("Error - Invalid slice region."); |
291 |
} |
292 |
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 |
} |
296 |
|
297 |
const DataTypes::ValueType& otherData=otherTemp->getVector(); |
298 |
const DataTypes::ShapeType& otherShape=otherTemp->getShape(); |
299 |
// copy slice from other default value to this default value |
300 |
// getDefaultValue().copySliceFrom(otherTemp->getDefaultValue(), regionLoopRange); |
301 |
DataTypes::copySliceFrom(m_data,getShape(),getDefaultOffset(),otherData,otherShape,otherTemp->getDefaultOffset(),regionLoopRange); |
302 |
|
303 |
// 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 |
addTag(pos->first); |
308 |
} |
309 |
} |
310 |
|
311 |
// loop through the tag values copying slices from other to this |
312 |
for (pos=m_offsetLookup.begin();pos!=m_offsetLookup.end();pos++) { |
313 |
// 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 |
} |
317 |
|
318 |
} |
319 |
|
320 |
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 |
if (dpno<0 || dpno>numDataPoints-1) { |
334 |
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 |
// void |
351 |
// DataTagged::setTaggedValues(const TagListType& tagKeys, |
352 |
// const ValueListType& values) |
353 |
// { |
354 |
// addTaggedValues(tagKeys,values); |
355 |
// } |
356 |
|
357 |
void |
358 |
DataTagged::setTaggedValue(int tagKey, |
359 |
const DataTypes::ShapeType& pointshape, |
360 |
const ValueType& value, |
361 |
int dataOffset) |
362 |
{ |
363 |
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 |
} |
379 |
|
380 |
|
381 |
/* |
382 |
void |
383 |
DataTagged::setTaggedValue(int tagKey, |
384 |
const DataArrayView& value) |
385 |
{ |
386 |
if (!getPointDataView().checkShape(value.getShape())) { |
387 |
throw DataException(DataTypes::createShapeErrorMessage( |
388 |
"Error - Cannot setTaggedValue due to shape mismatch.", value.getShape(),getShape())); |
389 |
} |
390 |
DataMapType::iterator pos(m_offsetLookup.find(tagKey)); |
391 |
if (pos==m_offsetLookup.end()) { |
392 |
// tag couldn't be found so use addTaggedValue |
393 |
addTaggedValue(tagKey,value); |
394 |
} else { |
395 |
// 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 |
} |
400 |
} |
401 |
}*/ |
402 |
|
403 |
// void |
404 |
// DataTagged::addTaggedValues(const TagListType& tagKeys, |
405 |
// const ValueListType& values) |
406 |
// { |
407 |
// if (values.size()==0) { |
408 |
// // copy the current default value for each of the tags |
409 |
// TagListType::const_iterator iT; |
410 |
// for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) { |
411 |
// // the point data view for DataTagged points at the default value |
412 |
// addTaggedValue(*iT,getPointDataView()); |
413 |
// } |
414 |
// } else if (values.size()==1 && tagKeys.size()>1) { |
415 |
// // assume the one given value will be used for all tag values |
416 |
// TagListType::const_iterator iT; |
417 |
// for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) { |
418 |
// addTaggedValue(*iT,values[0]); |
419 |
// } |
420 |
// } else { |
421 |
// if (tagKeys.size()!=values.size()) { |
422 |
// stringstream temp; |
423 |
// temp << "Error - (addTaggedValue) Number of tags: " << tagKeys.size() |
424 |
// << " doesn't match number of values: " << values.size(); |
425 |
// throw DataException(temp.str()); |
426 |
// } else { |
427 |
// unsigned int i; |
428 |
// for (i=0;i<tagKeys.size();i++) { |
429 |
// addTaggedValue(tagKeys[i],values[i]); |
430 |
// } |
431 |
// } |
432 |
// } |
433 |
// } |
434 |
|
435 |
|
436 |
void |
437 |
DataTagged::addTaggedValues(const TagListType& tagKeys, |
438 |
const ValueBatchType& values, |
439 |
const ShapeType& vShape) |
440 |
{ |
441 |
DataTypes::ValueType t(values.size(),0); |
442 |
for (size_t i=0;i<values.size();++i) |
443 |
{ |
444 |
t[i]=values[i]; |
445 |
} |
446 |
addTaggedValues(tagKeys,t,vShape); |
447 |
} |
448 |
|
449 |
|
450 |
// Note: The check to see if vShape==our shape is done in the addTaggedValue method |
451 |
void |
452 |
DataTagged::addTaggedValues(const TagListType& tagKeys, |
453 |
const ValueType& values, |
454 |
const ShapeType& vShape) |
455 |
{ |
456 |
int n=getNoValues(); |
457 |
int numVals=values.size()/n; |
458 |
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 |
addTag(*iT); |
464 |
} |
465 |
} else if (numVals==1 && tagKeys.size()>1) { |
466 |
// 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 |
addTaggedValue(*iT, vShape, values,0); |
470 |
} |
471 |
} else { |
472 |
if (tagKeys.size()!=numVals) { |
473 |
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 |
unsigned int i; |
479 |
int offset=0; |
480 |
for (i=0;i<tagKeys.size();i++ ,offset+=n) { |
481 |
addTaggedValue(tagKeys[i],vShape,values,offset); |
482 |
} |
483 |
} |
484 |
} |
485 |
} |
486 |
|
487 |
|
488 |
|
489 |
|
490 |
void |
491 |
DataTagged::addTaggedValue(int tagKey, |
492 |
const DataTypes::ShapeType& pointshape, |
493 |
const ValueType& value, |
494 |
int dataOffset) |
495 |
{ |
496 |
if (!DataTypes::checkShape(getShape(), pointshape)) { |
497 |
throw DataException(DataTypes::createShapeErrorMessage( |
498 |
"Error - Cannot addTaggedValue due to shape mismatch.", pointshape,getShape())); |
499 |
} |
500 |
DataMapType::iterator pos(m_offsetLookup.find(tagKey)); |
501 |
if (pos!=m_offsetLookup.end()) { |
502 |
// tag already exists so use setTaggedValue |
503 |
setTaggedValue(tagKey,pointshape, value, dataOffset); |
504 |
} 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 |
int newSize=m_data.size()+getNoValues(); |
513 |
m_data.resize(newSize,0.,newSize); |
514 |
for (int i=0;i<oldSize;i++) { |
515 |
m_data[i]=m_data_temp[i]; |
516 |
} |
517 |
for (int i=0;i<getNoValues();i++) { |
518 |
m_data[oldSize+i]=value[i+dataOffset]; |
519 |
} |
520 |
} |
521 |
} |
522 |
|
523 |
|
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 |
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 |
using namespace escript::DataTypes; |
602 |
string empty=""; |
603 |
stringstream temp; |
604 |
DataMapType::const_iterator i; |
605 |
temp << "Tag(Default)" << endl; |
606 |
temp << pointToString(m_data,getShape(),getDefaultOffset(),empty) << endl; |
607 |
// create a temporary view as the offset will be changed |
608 |
// DataArrayView tempView(getPointDataView().getData(), getPointDataView().getShape()); |
609 |
for (i=m_offsetLookup.begin();i!=m_offsetLookup.end();++i) { |
610 |
temp << "Tag(" << i->first << ")" << endl; |
611 |
temp << pointToString(m_data,getShape(),i->second,empty) << endl; |
612 |
// tempView.setOffset(i->second); |
613 |
// temp << tempView.toString() << endl; |
614 |
} |
615 |
return temp.str(); |
616 |
} |
617 |
|
618 |
DataTypes::ValueType::size_type |
619 |
DataTagged::getPointOffset(int sampleNo, |
620 |
int dataPointNo) const |
621 |
{ |
622 |
int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo); |
623 |
DataMapType::const_iterator pos(m_offsetLookup.find(tagKey)); |
624 |
DataTypes::ValueType::size_type offset=m_defaultValueOffset; |
625 |
if (pos!=m_offsetLookup.end()) { |
626 |
offset=pos->second; |
627 |
} |
628 |
return offset; |
629 |
} |
630 |
|
631 |
// 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 |
{ |
649 |
DataMapType::const_iterator pos(m_offsetLookup.find(tag)); |
650 |
DataTypes::ValueType::size_type offset=m_defaultValueOffset; |
651 |
if (pos!=m_offsetLookup.end()) { |
652 |
offset=pos->second; |
653 |
} |
654 |
return offset; |
655 |
} |
656 |
|
657 |
DataTypes::ValueType::const_reference |
658 |
DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i) const |
659 |
{ |
660 |
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 |
} |
670 |
|
671 |
|
672 |
DataTypes::ValueType::reference |
673 |
DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i) |
674 |
{ |
675 |
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 |
} |
685 |
|
686 |
|
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 |
void |
702 |
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 |
ValueType& evVec=temp_ev->getVector(); |
712 |
const ShapeType& evShape=temp_ev->getShape(); |
713 |
for (i=thisLookup.begin();i!=thisLookupEnd;i++) { |
714 |
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 |
} |
723 |
// symmetric(m_data,getShape(),getDefaultOffset(), |
724 |
DataMaths::symmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset()); |
725 |
} |
726 |
|
727 |
|
728 |
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 |
ValueType& evVec=temp_ev->getVector(); |
739 |
const ShapeType& evShape=temp_ev->getShape(); |
740 |
for (i=thisLookup.begin();i!=thisLookupEnd;i++) { |
741 |
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 |
} |
748 |
DataMaths::nonsymmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset()); |
749 |
} |
750 |
|
751 |
|
752 |
void |
753 |
DataTagged::trace(DataAbstract* ev, int axis_offset) |
754 |
{ |
755 |
DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev); |
756 |
if (temp_ev==0) { |
757 |
throw DataException("Error - DataTagged::trace casting to DataTagged failed (probably a programming error)."); |
758 |
} |
759 |
const DataTagged::DataMapType& thisLookup=getTagLookup(); |
760 |
DataTagged::DataMapType::const_iterator i; |
761 |
DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end(); |
762 |
ValueType& evVec=temp_ev->getVector(); |
763 |
const ShapeType& evShape=temp_ev->getShape(); |
764 |
for (i=thisLookup.begin();i!=thisLookupEnd;i++) { |
765 |
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 |
} |
772 |
DataMaths::trace(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset); |
773 |
} |
774 |
|
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 |
ValueType& evVec=temp_ev->getVector(); |
786 |
const ShapeType& evShape=temp_ev->getShape(); |
787 |
for (i=thisLookup.begin();i!=thisLookupEnd;i++) { |
788 |
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 |
} |
795 |
DataMaths::transpose(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset); |
796 |
} |
797 |
|
798 |
void |
799 |
DataTagged::swapaxes(DataAbstract* ev, int axis0, int axis1) |
800 |
{ |
801 |
DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev); |
802 |
if (temp_ev==0) { |
803 |
throw DataException("Error - DataTagged::swapaxes casting to DataTagged failed (probably a programming error)."); |
804 |
} |
805 |
const DataTagged::DataMapType& thisLookup=getTagLookup(); |
806 |
DataTagged::DataMapType::const_iterator i; |
807 |
DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end(); |
808 |
ValueType& evVec=temp_ev->getVector(); |
809 |
const ShapeType& evShape=temp_ev->getShape(); |
810 |
for (i=thisLookup.begin();i!=thisLookupEnd;i++) { |
811 |
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 |
} |
818 |
DataMaths::swapaxes(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis0,axis1); |
819 |
} |
820 |
|
821 |
void |
822 |
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 |
ValueType& evVec=temp_ev->getVector(); |
832 |
const ShapeType& evShape=temp_ev->getShape(); |
833 |
for (i=thisLookup.begin();i!=thisLookupEnd;i++) { |
834 |
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 |
} |
841 |
DataMaths::eigenvalues(m_data,getShape(),getDefaultOffset(),evVec, evShape, temp_ev->getDefaultOffset()); |
842 |
} |
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 |
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 |
for (i=thisLookup.begin();i!=thisLookupEnd;i++) { |
862 |
temp_ev->addTag(i->first); |
863 |
temp_V->addTag(i->first); |
864 |
/* DataArrayView thisView=getDataPointByTag(i->first); |
865 |
DataArrayView evView=temp_ev->getDataPointByTag(i->first); |
866 |
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 |
} |
874 |
DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),getDefaultOffset(),evVec, evShape, |
875 |
temp_ev->getDefaultOffset(),VVec,VShape, |
876 |
temp_V->getDefaultOffset(), tol); |
877 |
|
878 |
|
879 |
} |
880 |
|
881 |
void |
882 |
DataTagged::setToZero(){ |
883 |
DataTypes::ValueType::size_type n=m_data.size(); |
884 |
for (int i=0; i<n ;++i) m_data[i]=0.; |
885 |
} |
886 |
|
887 |
void |
888 |
DataTagged::dump(const std::string fileName) const |
889 |
{ |
890 |
#ifdef USE_NETCDF |
891 |
const int ldims=DataTypes::maxRank+1; |
892 |
const NcDim* ncdims[ldims]; |
893 |
NcVar *var, *tags_var; |
894 |
int rank = getRank(); |
895 |
int type= getFunctionSpace().getTypeCode(); |
896 |
int ndims =0; |
897 |
long dims[ldims]; |
898 |
const double* d_ptr=&(m_data[0]); |
899 |
DataTypes::ShapeType shape = getShape(); |
900 |
int mpi_iam=getFunctionSpace().getDomain()->getMPIRank(); |
901 |
int mpi_num=getFunctionSpace().getDomain()->getMPISize(); |
902 |
#ifdef PASO_MPI |
903 |
MPI_Status status; |
904 |
#endif |
905 |
|
906 |
#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 |
// netCDF error handler |
912 |
NcError err(NcError::verbose_nonfatal); |
913 |
// Create the file. |
914 |
char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam); |
915 |
NcFile dataFile(newFileName, NcFile::Replace); |
916 |
// check if writing was successful |
917 |
if (!dataFile.is_valid()) |
918 |
throw DataException("Error - DataTagged:: opening of netCDF file for output failed."); |
919 |
if (!dataFile.add_att("type_id",1) ) |
920 |
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 |
throw DataException("Error - DataTagged:: appending ncdimension 0 to netCDF file failed."); |
930 |
} |
931 |
if ( rank >1 ) { |
932 |
dims[1]=shape[1]; |
933 |
if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) ) |
934 |
throw DataException("Error - DataTagged:: appending ncdimension 1 to netCDF file failed."); |
935 |
} |
936 |
if ( rank >2 ) { |
937 |
dims[2]=shape[2]; |
938 |
if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) ) |
939 |
throw DataException("Error - DataTagged:: appending ncdimension 2 to netCDF file failed."); |
940 |
} |
941 |
if ( rank >3 ) { |
942 |
dims[3]=shape[3]; |
943 |
if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) ) |
944 |
throw DataException("Error - DataTagged:: appending ncdimension 3 to netCDF file failed."); |
945 |
} |
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 |
int* tags =(int*) esysUtils::malloc(ntags*sizeof(int)); |
952 |
int c=1; |
953 |
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 |
{ |
958 |
esysUtils::free(tags); |
959 |
throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed."); |
960 |
} |
961 |
if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) ) |
962 |
{ |
963 |
esysUtils::free(tags); |
964 |
throw DataException("Error - DataTagged:: appending tags to netCDF file failed."); |
965 |
} |
966 |
if (! (tags_var->put(tags,dims[rank])) ) |
967 |
{ |
968 |
esysUtils::free(tags); |
969 |
throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed."); |
970 |
} |
971 |
if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) ) |
972 |
{ |
973 |
esysUtils::free(tags); |
974 |
throw DataException("Error - DataTagged:: appending variable to netCDF file failed."); |
975 |
} |
976 |
if (! (var->put(d_ptr,dims)) ) |
977 |
{ |
978 |
esysUtils::free(tags); |
979 |
throw DataException("Error - DataTagged:: copy data to netCDF buffer failed."); |
980 |
} |
981 |
#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 |
#else |
985 |
throw DataException("Error - DataTagged:: dump is not configured with netCDF. Please contact your installation manager."); |
986 |
#endif |
987 |
} |
988 |
|
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 |
} // end of namespace |