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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2005 - (show annotations)
Mon Nov 10 01:21:39 2008 UTC (11 years ago) by jfenwick
File size: 35373 byte(s)
Bringing all changes across from schroedinger.
(Note this does not mean development is done, just that it will happen
on the trunk for now).
If anyone notices any problems please contact me.


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 : parent(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 : parent(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 : parent(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 : parent(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 : parent(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 : parent(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 : parent(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 unsigned int n=getNoValues();
457 unsigned 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 DataTypes::ValueType::size_type
632 DataTagged::getPointOffset(int sampleNo,
633 int dataPointNo)
634 {
635 int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
636 DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));
637 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
638 if (pos!=m_offsetLookup.end()) {
639 offset=pos->second;
640 }
641 return offset;
642 }
643
644 // DataArrayView
645 // DataTagged::getDataPointByTag(int tag) const
646 // {
647 // DataMapType::const_iterator pos(m_offsetLookup.find(tag));
648 // DataTypes::ValueType::size_type offset=m_defaultValueOffset;
649 // if (pos!=m_offsetLookup.end()) {
650 // offset=pos->second;
651 // }
652 // DataArrayView temp(getPointDataView());
653 // temp.setOffset(offset);
654 // return temp;
655 // }
656 //
657
658
659 DataTypes::ValueType::size_type
660 DataTagged::getOffsetForTag(int tag) const
661 {
662 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
663 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
664 if (pos!=m_offsetLookup.end()) {
665 offset=pos->second;
666 }
667 return offset;
668 }
669
670 DataTypes::ValueType::const_reference
671 DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i) const
672 {
673 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
674 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
675 if (pos!=m_offsetLookup.end()) {
676 offset=pos->second;
677 }
678 return m_data[offset+i];
679 /* DataArrayView temp(getPointDataView());
680 temp.setOffset(offset);
681 return temp.getData()[offset+i];*/
682 }
683
684
685 DataTypes::ValueType::reference
686 DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i)
687 {
688 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
689 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
690 if (pos!=m_offsetLookup.end()) {
691 offset=pos->second;
692 }
693 return m_data[offset+i];
694 /* DataArrayView temp(getPointDataView());
695 temp.setOffset(offset);
696 return temp.getData()[offset+i];*/
697 }
698
699
700
701
702
703
704 // DataArrayView
705 // DataTagged::getDataPoint(int sampleNo,
706 // int dataPointNo)
707 // {
708 // EsysAssert(validSampleNo(sampleNo),"(getDataPoint) Invalid sampleNo: " << sampleNo);
709 // int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
710 // return getDataPointByTag(tagKey);
711 // }
712
713
714 void
715 DataTagged::symmetric(DataAbstract* ev)
716 {
717 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
718 if (temp_ev==0) {
719 throw DataException("Error - DataTagged::symmetric casting to DataTagged failed (probably a programming error).");
720 }
721 const DataTagged::DataMapType& thisLookup=getTagLookup();
722 DataTagged::DataMapType::const_iterator i;
723 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
724 ValueType& evVec=temp_ev->getVector();
725 const ShapeType& evShape=temp_ev->getShape();
726 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
727 temp_ev->addTag(i->first);
728 DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
729 // DataArrayView thisView=getDataPointByTag(i->first);
730 // DataArrayView evView=temp_ev->getDataPointByTag(i->first);
731 DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
732
733 // DataArrayView::symmetric(thisView,0,evView,0);
734 DataMaths::symmetric(m_data,getShape(),offset,evVec, evShape, evoffset);
735 }
736 // symmetric(m_data,getShape(),getDefaultOffset(),
737 DataMaths::symmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
738 }
739
740
741 void
742 DataTagged::nonsymmetric(DataAbstract* ev)
743 {
744 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
745 if (temp_ev==0) {
746 throw DataException("Error - DataTagged::nonsymmetric casting to DataTagged failed (probably a programming error).");
747 }
748 const DataTagged::DataMapType& thisLookup=getTagLookup();
749 DataTagged::DataMapType::const_iterator i;
750 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
751 ValueType& evVec=temp_ev->getVector();
752 const ShapeType& evShape=temp_ev->getShape();
753 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
754 temp_ev->addTag(i->first);
755 /* DataArrayView thisView=getDataPointByTag(i->first);
756 DataArrayView evView=temp_ev->getDataPointByTag(i->first);*/
757 DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
758 DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
759 DataMaths::nonsymmetric(m_data,getShape(),offset,evVec, evShape, evoffset);
760 }
761 DataMaths::nonsymmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
762 }
763
764
765 void
766 DataTagged::trace(DataAbstract* ev, int axis_offset)
767 {
768 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
769 if (temp_ev==0) {
770 throw DataException("Error - DataTagged::trace casting to DataTagged failed (probably a programming error).");
771 }
772 const DataTagged::DataMapType& thisLookup=getTagLookup();
773 DataTagged::DataMapType::const_iterator i;
774 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
775 ValueType& evVec=temp_ev->getVector();
776 const ShapeType& evShape=temp_ev->getShape();
777 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
778 temp_ev->addTag(i->first);
779 // DataArrayView thisView=getDataPointByTag(i->first);
780 // DataArrayView evView=temp_ev->getDataPointByTag(i->first);
781 DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
782 DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
783 DataMaths::trace(m_data,getShape(),offset,evVec, evShape, evoffset, axis_offset);
784 }
785 DataMaths::trace(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
786 }
787
788 void
789 DataTagged::transpose(DataAbstract* ev, int axis_offset)
790 {
791 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
792 if (temp_ev==0) {
793 throw DataException("Error - DataTagged::transpose casting to DataTagged failed (probably a programming error).");
794 }
795 const DataTagged::DataMapType& thisLookup=getTagLookup();
796 DataTagged::DataMapType::const_iterator i;
797 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
798 ValueType& evVec=temp_ev->getVector();
799 const ShapeType& evShape=temp_ev->getShape();
800 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
801 temp_ev->addTag(i->first);
802 // DataArrayView thisView=getDataPointByTag(i->first);
803 // DataArrayView evView=temp_ev->getDataPointByTag(i->first);
804 DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
805 DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
806 DataMaths::transpose(m_data,getShape(),offset,evVec, evShape, evoffset, axis_offset);
807 }
808 DataMaths::transpose(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
809 }
810
811 void
812 DataTagged::swapaxes(DataAbstract* ev, int axis0, int axis1)
813 {
814 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
815 if (temp_ev==0) {
816 throw DataException("Error - DataTagged::swapaxes casting to DataTagged failed (probably a programming error).");
817 }
818 const DataTagged::DataMapType& thisLookup=getTagLookup();
819 DataTagged::DataMapType::const_iterator i;
820 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
821 ValueType& evVec=temp_ev->getVector();
822 const ShapeType& evShape=temp_ev->getShape();
823 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
824 temp_ev->addTag(i->first);
825 /* DataArrayView thisView=getDataPointByTag(i->first);
826 DataArrayView evView=temp_ev->getDataPointByTag(i->first);*/
827 DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
828 DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
829 DataMaths::swapaxes(m_data,getShape(),offset,evVec, evShape, evoffset,axis0,axis1);
830 }
831 DataMaths::swapaxes(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis0,axis1);
832 }
833
834 void
835 DataTagged::eigenvalues(DataAbstract* ev)
836 {
837 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
838 if (temp_ev==0) {
839 throw DataException("Error - DataTagged::eigenvalues casting to DataTagged failed (propably a programming error).");
840 }
841 const DataTagged::DataMapType& thisLookup=getTagLookup();
842 DataTagged::DataMapType::const_iterator i;
843 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
844 ValueType& evVec=temp_ev->getVector();
845 const ShapeType& evShape=temp_ev->getShape();
846 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
847 temp_ev->addTag(i->first);
848 // DataArrayView thisView=getDataPointByTag(i->first);
849 // DataArrayView evView=temp_ev->getDataPointByTag(i->first);
850 DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
851 DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
852 DataMaths::eigenvalues(m_data,getShape(),offset,evVec, evShape, evoffset);
853 }
854 DataMaths::eigenvalues(m_data,getShape(),getDefaultOffset(),evVec, evShape, temp_ev->getDefaultOffset());
855 }
856 void
857 DataTagged::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
858 {
859 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
860 if (temp_ev==0) {
861 throw DataException("Error - DataTagged::eigenvalues_and_eigenvectors casting to DataTagged failed (propably a programming error).");
862 }
863 DataTagged* temp_V=dynamic_cast<DataTagged*>(V);
864 if (temp_V==0) {
865 throw DataException("Error - DataTagged::eigenvalues_and_eigenvectors casting to DataTagged failed (propably a programming error).");
866 }
867 const DataTagged::DataMapType& thisLookup=getTagLookup();
868 DataTagged::DataMapType::const_iterator i;
869 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
870 ValueType& evVec=temp_ev->getVector();
871 const ShapeType& evShape=temp_ev->getShape();
872 ValueType& VVec=temp_V->getVector();
873 const ShapeType& VShape=temp_V->getShape();
874 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
875 temp_ev->addTag(i->first);
876 temp_V->addTag(i->first);
877 /* DataArrayView thisView=getDataPointByTag(i->first);
878 DataArrayView evView=temp_ev->getDataPointByTag(i->first);
879 DataArrayView VView=temp_V->getDataPointByTag(i->first);*/
880 DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
881 DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
882 DataTypes::ValueType::size_type Voffset=temp_V->getOffsetForTag(i->first);
883 /* DataArrayView::eigenvalues_and_eigenvectors(thisView,0,evView,0,VView,0,tol);*/
884 DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),offset,evVec, evShape, evoffset,VVec,VShape,Voffset,tol);
885
886 }
887 DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),getDefaultOffset(),evVec, evShape,
888 temp_ev->getDefaultOffset(),VVec,VShape,
889 temp_V->getDefaultOffset(), tol);
890
891
892 }
893
894 void
895 DataTagged::setToZero(){
896 DataTypes::ValueType::size_type n=m_data.size();
897 for (int i=0; i<n ;++i) m_data[i]=0.;
898 }
899
900 void
901 DataTagged::dump(const std::string fileName) const
902 {
903 #ifdef USE_NETCDF
904 const int ldims=DataTypes::maxRank+1;
905 const NcDim* ncdims[ldims];
906 NcVar *var, *tags_var;
907 int rank = getRank();
908 int type= getFunctionSpace().getTypeCode();
909 int ndims =0;
910 long dims[ldims];
911 const double* d_ptr=&(m_data[0]);
912 DataTypes::ShapeType shape = getShape();
913 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
914 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
915 #ifdef PASO_MPI
916 MPI_Status status;
917 #endif
918
919 #ifdef PASO_MPI
920 /* Serialize NetCDF I/O */
921 if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81803, MPI_COMM_WORLD, &status);
922 #endif
923
924 // netCDF error handler
925 NcError err(NcError::verbose_nonfatal);
926 // Create the file.
927 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
928 NcFile dataFile(newFileName, NcFile::Replace);
929 // check if writing was successful
930 if (!dataFile.is_valid())
931 throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");
932 if (!dataFile.add_att("type_id",1) )
933 throw DataException("Error - DataTagged:: appending data type to netCDF file failed.");
934 if (!dataFile.add_att("rank",rank) )
935 throw DataException("Error - DataTagged:: appending rank attribute to netCDF file failed.");
936 if (!dataFile.add_att("function_space_type",type))
937 throw DataException("Error - DataTagged:: appending function space attribute to netCDF file failed.");
938 ndims=rank+1;
939 if ( rank >0 ) {
940 dims[0]=shape[0];
941 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
942 throw DataException("Error - DataTagged:: appending ncdimension 0 to netCDF file failed.");
943 }
944 if ( rank >1 ) {
945 dims[1]=shape[1];
946 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
947 throw DataException("Error - DataTagged:: appending ncdimension 1 to netCDF file failed.");
948 }
949 if ( rank >2 ) {
950 dims[2]=shape[2];
951 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
952 throw DataException("Error - DataTagged:: appending ncdimension 2 to netCDF file failed.");
953 }
954 if ( rank >3 ) {
955 dims[3]=shape[3];
956 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
957 throw DataException("Error - DataTagged:: appending ncdimension 3 to netCDF file failed.");
958 }
959 const DataTagged::DataMapType& thisLookup=getTagLookup();
960 DataTagged::DataMapType::const_iterator i;
961 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
962 int ntags=1;
963 for (i=thisLookup.begin();i!=thisLookupEnd;i++) ntags++;
964 int* tags =(int*) esysUtils::malloc(ntags*sizeof(int));
965 int c=1;
966 tags[0]=-1;
967 for (i=thisLookup.begin();i!=thisLookupEnd;i++) tags[c++]=i->first;
968 dims[rank]=ntags;
969 if (! (ncdims[rank] = dataFile.add_dim("num_tags", dims[rank])) )
970 {
971 esysUtils::free(tags);
972 throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");
973 }
974 if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )
975 {
976 esysUtils::free(tags);
977 throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");
978 }
979 if (! (tags_var->put(tags,dims[rank])) )
980 {
981 esysUtils::free(tags);
982 throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");
983 }
984 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
985 {
986 esysUtils::free(tags);
987 throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");
988 }
989 if (! (var->put(d_ptr,dims)) )
990 {
991 esysUtils::free(tags);
992 throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");
993 }
994 #ifdef PASO_MPI
995 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81803, MPI_COMM_WORLD);
996 #endif
997 #else
998 throw DataException("Error - DataTagged:: dump is not configured with netCDF. Please contact your installation manager.");
999 #endif
1000 }
1001
1002 DataTypes::ValueType&
1003 DataTagged::getVector()
1004 {
1005 return m_data;
1006 }
1007
1008 const DataTypes::ValueType&
1009 DataTagged::getVector() const
1010 {
1011 return m_data;
1012 }
1013
1014 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26