/[escript]/branches/arrexp_2137_win/escript/src/DataTagged.cpp
ViewVC logotype

Contents of /branches/arrexp_2137_win/escript/src/DataTagged.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2202 - (show annotations)
Fri Jan 9 01:28:32 2009 UTC (10 years, 5 months ago) by jfenwick
File size: 35400 byte(s)
Branching the array experiments from version 2137.
The idea is to make the changes required for the c++ changes to compile 
on windows without bringing in the later python changes.


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 (unsigned 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 (unsigned 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 (unsigned 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