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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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_data(other.m_data),
150 m_offsetLookup(other.m_offsetLookup)
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 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 DataAbstract*
227 DataTagged::deepCopy()
228 {
229 return new DataTagged(*this);
230 }
231
232 DataAbstract*
233 DataTagged::getSlice(const DataTypes::RegionType& region) const
234 {
235 return new DataTagged(*this, region);
236 }
237
238 DataTagged::DataTagged(const DataTagged& other,
239 const DataTypes::RegionType& region)
240 : DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
241 {
242 // slice constructor
243
244 // get the shape of the slice to copy from other
245 DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
246 DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
247
248 // allocate enough space in this for all values
249 // (need to add one to allow for the default value)
250 int len = DataTypes::noValues(regionShape)*(other.m_offsetLookup.size()+1);
251 m_data.resize(len,0.0,len);
252
253 // create the data view
254 // DataArrayView temp(m_data,regionShape);
255 // setPointDataView(temp);
256
257 // copy the default value from other to this
258 // 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
263 // loop through the tag values copying these
264 DataMapType::const_iterator pos;
265 DataTypes::ValueType::size_type tagOffset=getNoValues();
266 for (pos=other.m_offsetLookup.begin();pos!=other.m_offsetLookup.end();pos++){
267 // getPointDataView().copySlice(tagOffset,other.getPointDataView(),pos->second,regionLoopRange);
268 DataTypes::copySlice(m_data,getShape(),tagOffset,otherData, otherShape, pos->second, regionLoopRange);
269 m_offsetLookup.insert(DataMapType::value_type(pos->first,tagOffset));
270 tagOffset+=getNoValues();
271 }
272 }
273
274 void
275 DataTagged::setSlice(const DataAbstract* other,
276 const DataTypes::RegionType& region)
277 {
278
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 throw DataException("Programming error - casting to DataTagged.");
284 }
285
286 // determine shape of the specified region
287 DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
288
289 // modify region specification as needed to match rank of this object
290 DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
291
292 // ensure rank/shape of this object is compatible with specified region
293 if (getRank()!=region.size()) {
294 throw DataException("Error - Invalid slice region.");
295 }
296 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 }
300
301 const DataTypes::ValueType& otherData=otherTemp->getVector();
302 const DataTypes::ShapeType& otherShape=otherTemp->getShape();
303 // copy slice from other default value to this default value
304 // getDefaultValue().copySliceFrom(otherTemp->getDefaultValue(), regionLoopRange);
305 DataTypes::copySliceFrom(m_data,getShape(),getDefaultOffset(),otherData,otherShape,otherTemp->getDefaultOffset(),regionLoopRange);
306
307 // 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 addTag(pos->first);
312 }
313 }
314
315 // loop through the tag values copying slices from other to this
316 for (pos=m_offsetLookup.begin();pos!=m_offsetLookup.end();pos++) {
317 // 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 }
321
322 }
323
324 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 if (dpno<0 || dpno>numDataPoints-1) {
338 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 // void
355 // DataTagged::setTaggedValues(const TagListType& tagKeys,
356 // const ValueListType& values)
357 // {
358 // addTaggedValues(tagKeys,values);
359 // }
360
361 void
362 DataTagged::setTaggedValue(int tagKey,
363 const DataTypes::ShapeType& pointshape,
364 const ValueType& value,
365 int dataOffset)
366 {
367 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 }
383
384
385 /*
386 void
387 DataTagged::setTaggedValue(int tagKey,
388 const DataArrayView& value)
389 {
390 if (!getPointDataView().checkShape(value.getShape())) {
391 throw DataException(DataTypes::createShapeErrorMessage(
392 "Error - Cannot setTaggedValue due to shape mismatch.", value.getShape(),getShape()));
393 }
394 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
395 if (pos==m_offsetLookup.end()) {
396 // tag couldn't be found so use addTaggedValue
397 addTaggedValue(tagKey,value);
398 } else {
399 // 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 }
404 }
405 }*/
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 for (size_t i=0;i<values.size();++i)
447 {
448 t[i]=values[i];
449 }
450 addTaggedValues(tagKeys,t,vShape);
451 }
452
453
454 // Note: The check to see if vShape==our shape is done in the addTaggedValue method
455 void
456 DataTagged::addTaggedValues(const TagListType& tagKeys,
457 const ValueType& values,
458 const ShapeType& vShape)
459 {
460 int n=getNoValues();
461 int numVals=values.size()/n;
462 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 addTag(*iT);
468 }
469 } else if (numVals==1 && tagKeys.size()>1) {
470 // 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 addTaggedValue(*iT, vShape, values,0);
474 }
475 } else {
476 if (tagKeys.size()!=numVals) {
477 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 unsigned int i;
483 int offset=0;
484 for (i=0;i<tagKeys.size();i++ ,offset+=n) {
485 addTaggedValue(tagKeys[i],vShape,values,offset);
486 }
487 }
488 }
489 }
490
491
492
493
494 void
495 DataTagged::addTaggedValue(int tagKey,
496 const DataTypes::ShapeType& pointshape,
497 const ValueType& value,
498 int dataOffset)
499 {
500 if (!DataTypes::checkShape(getShape(), pointshape)) {
501 throw DataException(DataTypes::createShapeErrorMessage(
502 "Error - Cannot addTaggedValue due to shape mismatch.", pointshape,getShape()));
503 }
504 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
505 if (pos!=m_offsetLookup.end()) {
506 // tag already exists so use setTaggedValue
507 setTaggedValue(tagKey,pointshape, value, dataOffset);
508 } 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 int newSize=m_data.size()+getNoValues();
517 m_data.resize(newSize,0.,newSize);
518 for (int i=0;i<oldSize;i++) {
519 m_data[i]=m_data_temp[i];
520 }
521 for (int i=0;i<getNoValues();i++) {
522 m_data[oldSize+i]=value[i+dataOffset];
523 }
524 }
525 }
526
527
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 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 using namespace escript::DataTypes;
606 string empty="";
607 stringstream temp;
608 DataMapType::const_iterator i;
609 temp << "Tag(Default)" << endl;
610 temp << pointToString(m_data,getShape(),getDefaultOffset(),empty) << endl;
611 // create a temporary view as the offset will be changed
612 // DataArrayView tempView(getPointDataView().getData(), getPointDataView().getShape());
613 for (i=m_offsetLookup.begin();i!=m_offsetLookup.end();++i) {
614 temp << "Tag(" << i->first << ")" << endl;
615 temp << pointToString(m_data,getShape(),i->second,empty) << endl;
616 // tempView.setOffset(i->second);
617 // temp << tempView.toString() << endl;
618 }
619 return temp.str();
620 }
621
622 DataTypes::ValueType::size_type
623 DataTagged::getPointOffset(int sampleNo,
624 int dataPointNo) const
625 {
626 int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
627 DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));
628 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
629 if (pos!=m_offsetLookup.end()) {
630 offset=pos->second;
631 }
632 return offset;
633 }
634
635 // 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 {
653 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
654 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
655 if (pos!=m_offsetLookup.end()) {
656 offset=pos->second;
657 }
658 return offset;
659 }
660
661 DataTypes::ValueType::const_reference
662 DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i) const
663 {
664 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 }
674
675
676 DataTypes::ValueType::reference
677 DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i)
678 {
679 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 }
689
690
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 void
706 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 ValueType& evVec=temp_ev->getVector();
716 const ShapeType& evShape=temp_ev->getShape();
717 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
718 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 }
727 // symmetric(m_data,getShape(),getDefaultOffset(),
728 DataMaths::symmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
729 }
730
731
732 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 ValueType& evVec=temp_ev->getVector();
743 const ShapeType& evShape=temp_ev->getShape();
744 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
745 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 }
752 DataMaths::nonsymmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
753 }
754
755
756 void
757 DataTagged::trace(DataAbstract* ev, int axis_offset)
758 {
759 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
760 if (temp_ev==0) {
761 throw DataException("Error - DataTagged::trace casting to DataTagged failed (probably a programming error).");
762 }
763 const DataTagged::DataMapType& thisLookup=getTagLookup();
764 DataTagged::DataMapType::const_iterator i;
765 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
766 ValueType& evVec=temp_ev->getVector();
767 const ShapeType& evShape=temp_ev->getShape();
768 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
769 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 }
776 DataMaths::trace(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
777 }
778
779 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 ValueType& evVec=temp_ev->getVector();
790 const ShapeType& evShape=temp_ev->getShape();
791 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
792 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 }
799 DataMaths::transpose(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
800 }
801
802 void
803 DataTagged::swapaxes(DataAbstract* ev, int axis0, int axis1)
804 {
805 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
806 if (temp_ev==0) {
807 throw DataException("Error - DataTagged::swapaxes casting to DataTagged failed (probably a programming error).");
808 }
809 const DataTagged::DataMapType& thisLookup=getTagLookup();
810 DataTagged::DataMapType::const_iterator i;
811 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
812 ValueType& evVec=temp_ev->getVector();
813 const ShapeType& evShape=temp_ev->getShape();
814 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
815 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 }
822 DataMaths::swapaxes(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis0,axis1);
823 }
824
825 void
826 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 ValueType& evVec=temp_ev->getVector();
836 const ShapeType& evShape=temp_ev->getShape();
837 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
838 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 }
845 DataMaths::eigenvalues(m_data,getShape(),getDefaultOffset(),evVec, evShape, temp_ev->getDefaultOffset());
846 }
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 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 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
866 temp_ev->addTag(i->first);
867 temp_V->addTag(i->first);
868 /* DataArrayView thisView=getDataPointByTag(i->first);
869 DataArrayView evView=temp_ev->getDataPointByTag(i->first);
870 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 }
878 DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),getDefaultOffset(),evVec, evShape,
879 temp_ev->getDefaultOffset(),VVec,VShape,
880 temp_V->getDefaultOffset(), tol);
881
882
883 }
884
885 void
886 DataTagged::setToZero(){
887 DataTypes::ValueType::size_type n=m_data.size();
888 for (int i=0; i<n ;++i) m_data[i]=0.;
889 }
890
891 void
892 DataTagged::dump(const std::string fileName) const
893 {
894 #ifdef USE_NETCDF
895 const int ldims=DataTypes::maxRank+1;
896 const NcDim* ncdims[ldims];
897 NcVar *var, *tags_var;
898 int rank = getRank();
899 int type= getFunctionSpace().getTypeCode();
900 int ndims =0;
901 long dims[ldims];
902 const double* d_ptr=&(m_data[0]);
903 DataTypes::ShapeType shape = getShape();
904 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
905 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
906 #ifdef PASO_MPI
907 MPI_Status status;
908 #endif
909
910 #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 // netCDF error handler
916 NcError err(NcError::verbose_nonfatal);
917 // Create the file.
918 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
919 NcFile dataFile(newFileName, NcFile::Replace);
920 // check if writing was successful
921 if (!dataFile.is_valid())
922 throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");
923 if (!dataFile.add_att("type_id",1) )
924 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 throw DataException("Error - DataTagged:: appending ncdimension 0 to netCDF file failed.");
934 }
935 if ( rank >1 ) {
936 dims[1]=shape[1];
937 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
938 throw DataException("Error - DataTagged:: appending ncdimension 1 to netCDF file failed.");
939 }
940 if ( rank >2 ) {
941 dims[2]=shape[2];
942 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
943 throw DataException("Error - DataTagged:: appending ncdimension 2 to netCDF file failed.");
944 }
945 if ( rank >3 ) {
946 dims[3]=shape[3];
947 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
948 throw DataException("Error - DataTagged:: appending ncdimension 3 to netCDF file failed.");
949 }
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 int* tags =(int*) esysUtils::malloc(ntags*sizeof(int));
956 int c=1;
957 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 {
962 esysUtils::free(tags);
963 throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");
964 }
965 if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )
966 {
967 esysUtils::free(tags);
968 throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");
969 }
970 if (! (tags_var->put(tags,dims[rank])) )
971 {
972 esysUtils::free(tags);
973 throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");
974 }
975 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
976 {
977 esysUtils::free(tags);
978 throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");
979 }
980 if (! (var->put(d_ptr,dims)) )
981 {
982 esysUtils::free(tags);
983 throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");
984 }
985 #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 #else
989 throw DataException("Error - DataTagged:: dump is not configured with netCDF. Please contact your installation manager.");
990 #endif
991 }
992
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 } // end of namespace

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26