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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1710 - (show annotations)
Fri Aug 15 06:24:20 2008 UTC (12 years, 6 months ago) by jfenwick
File size: 24611 byte(s)
Branch commit.

Bug fix


1
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 #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
25 using namespace std;
26
27 namespace escript {
28
29 DataTagged::DataTagged()
30 : DataAbstract(FunctionSpace())
31 {
32 // default constructor
33
34 // create a scalar default value
35 m_data.resize(1,0.,1);
36 DataArrayView temp(m_data,DataTypes::ShapeType());
37 setPointDataView(temp);
38 }
39
40 DataTagged::DataTagged(const TagListType& tagKeys,
41 const ValueListType& values,
42 const DataArrayView& defaultValue,
43 const FunctionSpace& what)
44 : DataAbstract(what)
45 {
46 // constructor
47
48 // initialise the array of data values
49 // the default value is always the first item in the values list
50 int len = defaultValue.noValues();
51 m_data.resize(len,0.,len);
52 for (int i=0; i<defaultValue.noValues(); i++) {
53 m_data[i]=defaultValue.getData(i);
54 }
55
56 // create the data view
57 DataArrayView temp(m_data,defaultValue.getShape());
58 setPointDataView(temp);
59
60 // add remaining tags and values
61 addTaggedValues(tagKeys,values);
62 }
63
64 DataTagged::DataTagged(const FunctionSpace& what,
65 const DataTypes::ShapeType &shape,
66 const int tags[],
67 const ValueType& data)
68 : DataAbstract(what)
69 {
70 // alternative constructor
71 // not unit_tested tested yet
72
73 // copy the data
74 m_data=data;
75
76 // create the view of the data
77 DataArrayView tempView(m_data,shape);
78 setPointDataView(tempView);
79
80 // create the tag lookup map
81 for (int sampleNo=0; sampleNo<getNumSamples(); sampleNo++) {
82 m_offsetLookup.insert(DataMapType::value_type(sampleNo,tags[sampleNo]));
83 }
84 }
85
86 DataTagged::DataTagged(const FunctionSpace& what,
87 const DataTypes::ShapeType &shape,
88 const TagListType& tags,
89 const ValueType& data)
90 : DataAbstract(what)
91 {
92 // alternative constructor
93 // not unit_tested tested yet
94
95 // copy the data
96 m_data=data;
97
98 // create the view of the data
99 DataArrayView tempView(m_data,shape);
100 setPointDataView(tempView);
101
102 // create the tag lookup map
103 for (int sampleNo=0; sampleNo<getNumSamples(); sampleNo++) {
104 m_offsetLookup.insert(DataMapType::value_type(sampleNo,tags[sampleNo]));
105 }
106 }
107
108
109 DataTagged::DataTagged(const DataTagged& other)
110 : DataAbstract(other.getFunctionSpace()),
111 m_data(other.m_data),
112 m_offsetLookup(other.m_offsetLookup)
113 {
114 // copy constructor
115
116 // create the data view
117 DataArrayView temp(m_data,other.getPointDataView().getShape());
118 setPointDataView(temp);
119 }
120
121 DataTagged::DataTagged(const DataConstant& other)
122 : DataAbstract(other.getFunctionSpace())
123 {
124 // copy constructor
125
126 // fill the default value with the constant value item from "other"
127 const DataArrayView& value=other.getPointDataView();
128 int len = value.noValues();
129 m_data.resize(len,0.,len);
130 for (int i=0; i<value.noValues(); i++) {
131 m_data[i]=value.getData(i);
132 }
133
134 // create the data view
135 DataArrayView temp(m_data,value.getShape());
136 setPointDataView(temp);
137 }
138
139 DataAbstract*
140 DataTagged::getSlice(const DataTypes::RegionType& region) const
141 {
142 return new DataTagged(*this, region);
143 }
144
145 DataTagged::DataTagged(const DataTagged& other,
146 const DataTypes::RegionType& region)
147 : DataAbstract(other.getFunctionSpace())
148 {
149 // slice constructor
150
151 // get the shape of the slice to copy from other
152 DataTypes::ShapeType regionShape(DataArrayView::getResultSliceShape(region));
153 DataTypes::RegionLoopRangeType regionLoopRange=getSliceRegionLoopRange(region);
154
155 // allocate enough space in this for all values
156 // (need to add one to allow for the default value)
157 int len = DataTypes::noValues(regionShape)*(other.m_offsetLookup.size()+1);
158 m_data.resize(len,0.0,len);
159
160 // create the data view
161 DataArrayView temp(m_data,regionShape);
162 setPointDataView(temp);
163
164 // copy the default value from other to this
165 getDefaultValue().copySlice(other.getDefaultValue(), regionLoopRange);
166
167 // loop through the tag values copying these
168 DataMapType::const_iterator pos;
169 DataTypes::ValueType::size_type tagOffset=getPointDataView().noValues();
170 for (pos=other.m_offsetLookup.begin();pos!=other.m_offsetLookup.end();pos++){
171 getPointDataView().copySlice(tagOffset,other.getPointDataView(),pos->second,regionLoopRange);
172 m_offsetLookup.insert(DataMapType::value_type(pos->first,tagOffset));
173 tagOffset+=getPointDataView().noValues();
174 }
175 }
176
177 void
178 DataTagged::setSlice(const DataAbstract* other,
179 const DataTypes::RegionType& region)
180 {
181
182 // other must be another DataTagged object
183 // Data:setSlice implementation should ensure this
184 const DataTagged* otherTemp=dynamic_cast<const DataTagged*>(other);
185 if (otherTemp==0) {
186 throw DataException("Programming error - casting to DataTagged.");
187 }
188
189 // determine shape of the specified region
190 DataTypes::ShapeType regionShape(DataArrayView::getResultSliceShape(region));
191
192 // modify region specification as needed to match rank of this object
193 DataTypes::RegionLoopRangeType regionLoopRange=getSliceRegionLoopRange(region);
194
195 // ensure rank/shape of this object is compatible with specified region
196 if (getPointDataView().getRank()!=region.size()) {
197 throw DataException("Error - Invalid slice region.");
198 }
199 if (otherTemp->getPointDataView().getRank()>0 && !other->getPointDataView().checkShape(regionShape)) {
200 throw DataException (other->getPointDataView().createShapeErrorMessage(
201 "Error - Couldn't copy slice due to shape mismatch.",regionShape));
202 }
203
204 // copy slice from other default value to this default value
205 getDefaultValue().copySliceFrom(otherTemp->getDefaultValue(), regionLoopRange);
206
207 // loop through tag values in other, adding any which aren't in this, using default value
208 DataMapType::const_iterator pos;
209 for (pos=otherTemp->m_offsetLookup.begin();pos!=otherTemp->m_offsetLookup.end();pos++) {
210 if (!isCurrentTag(pos->first)) {
211 addTag(pos->first);
212 }
213 }
214
215 // loop through the tag values copying slices from other to this
216 for (pos=m_offsetLookup.begin();pos!=m_offsetLookup.end();pos++) {
217 getDataPointByTag(pos->first).copySliceFrom(otherTemp->getDataPointByTag(pos->first), regionLoopRange);
218 }
219
220 }
221
222 int
223 DataTagged::getTagNumber(int dpno)
224 {
225 //
226 // Get the number of samples and data-points per sample
227 int numSamples = getNumSamples();
228 int numDataPointsPerSample = getNumDPPSample();
229 int numDataPoints = numSamples * numDataPointsPerSample;
230
231 if (numDataPointsPerSample==0) {
232 throw DataException("DataTagged::getTagNumber error: no data-points associated with this object.");
233 }
234
235 if (dpno<0 || dpno>numDataPoints-1) {
236 throw DataException("DataTagged::getTagNumber error: invalid data-point number supplied.");
237 }
238
239 //
240 // Determine the sample number which corresponds to this data-point number
241 int sampleNo = dpno / numDataPointsPerSample;
242
243 //
244 // Determine the tag number which corresponds to this sample number
245 int tagNo = getFunctionSpace().getTagFromSampleNo(sampleNo);
246
247 //
248 // return the tag number
249 return(tagNo);
250 }
251
252 void
253 DataTagged::setTaggedValues(const TagListType& tagKeys,
254 const ValueListType& values)
255 {
256 addTaggedValues(tagKeys,values);
257 }
258
259 void
260 DataTagged::setTaggedValue(int tagKey,
261 const DataArrayView& value)
262 {
263 if (!getPointDataView().checkShape(value.getShape())) {
264 throw DataException(getPointDataView().createShapeErrorMessage(
265 "Error - Cannot setTaggedValue due to shape mismatch.", value.getShape()));
266 }
267 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
268 if (pos==m_offsetLookup.end()) {
269 // tag couldn't be found so use addTaggedValue
270 addTaggedValue(tagKey,value);
271 } else {
272 // copy the values into the data array at the offset determined by m_offsetLookup
273 int offset=pos->second;
274 for (int i=0; i<getPointDataView().noValues(); i++) {
275 m_data[offset+i]=value.getData(i);
276 }
277 }
278 }
279
280 void
281 DataTagged::addTaggedValues(const TagListType& tagKeys,
282 const ValueListType& values)
283 {
284 if (values.size()==0) {
285 // copy the current default value for each of the tags
286 TagListType::const_iterator iT;
287 for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
288 // the point data view for DataTagged points at the default value
289 addTaggedValue(*iT,getPointDataView());
290 }
291 } else if (values.size()==1 && tagKeys.size()>1) {
292 // assume the one given value will be used for all tag values
293 TagListType::const_iterator iT;
294 for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
295 addTaggedValue(*iT,values[0]);
296 }
297 } else {
298 if (tagKeys.size()!=values.size()) {
299 stringstream temp;
300 temp << "Error - (addTaggedValue) Number of tags: " << tagKeys.size()
301 << " doesn't match number of values: " << values.size();
302 throw DataException(temp.str());
303 } else {
304 unsigned int i;
305 for (i=0;i<tagKeys.size();i++) {
306 addTaggedValue(tagKeys[i],values[i]);
307 }
308 }
309 }
310 }
311
312 void
313 DataTagged::addTaggedValue(int tagKey,
314 const DataArrayView& value)
315 {
316 if (!getPointDataView().checkShape(value.getShape())) {
317 throw DataException(getPointDataView().createShapeErrorMessage(
318 "Error - Cannot addTaggedValue due to shape mismatch.", value.getShape()));
319 }
320 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
321 if (pos!=m_offsetLookup.end()) {
322 // tag already exists so use setTaggedValue
323 setTaggedValue(tagKey,value);
324 } else {
325 // save the key and the location of its data in the lookup tab
326 m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));
327 // add the data given in "value" at the end of m_data
328 // need to make a temp copy of m_data, resize m_data, then copy
329 // all the old values plus the value to be added back into m_data
330 ValueType m_data_temp(m_data);
331 int oldSize=m_data.size();
332 int newSize=m_data.size()+value.noValues();
333 m_data.resize(newSize,0.,newSize);
334 for (int i=0;i<oldSize;i++) {
335 m_data[i]=m_data_temp[i];
336 }
337 for (int i=0;i<value.noValues();i++) {
338 m_data[oldSize+i]=value.getData(i);
339 }
340 }
341 }
342
343
344 void
345 DataTagged::addTag(int tagKey)
346 {
347 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
348 if (pos!=m_offsetLookup.end()) {
349 // tag already exists so use setTaggedValue
350 // setTaggedValue(tagKey,value);
351 } else {
352 // save the key and the location of its data in the lookup tab
353 m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));
354 // add the data given in "value" at the end of m_data
355 // need to make a temp copy of m_data, resize m_data, then copy
356 // all the old values plus the value to be added back into m_data
357 ValueType m_data_temp(m_data);
358 int oldSize=m_data.size();
359 int newSize=m_data.size()+getNoValues();
360 m_data.resize(newSize,0.,newSize);
361 for (int i=0;i<oldSize;i++) {
362 m_data[i]=m_data_temp[i];
363 }
364 for (int i=0;i<getNoValues();i++) {
365 m_data[oldSize+i]=m_data[m_defaultValueOffset+i];
366 }
367 }
368 }
369
370
371 double*
372 DataTagged::getSampleDataByTag(int tag)
373 {
374 DataMapType::iterator pos(m_offsetLookup.find(tag));
375 if (pos==m_offsetLookup.end()) {
376 // tag couldn't be found so return the default value
377 return &(m_data[0]);
378 } else {
379 // return the data-point corresponding to the given tag
380 return &(m_data[pos->second]);
381 }
382 }
383
384 string
385 DataTagged::toString() const
386 {
387 stringstream temp;
388 DataMapType::const_iterator i;
389 temp << "Tag(Default)" << endl;
390 temp << getDefaultValue().toString() << endl;
391 // create a temporary view as the offset will be changed
392 DataArrayView tempView(getPointDataView().getData(), getPointDataView().getShape());
393 for (i=m_offsetLookup.begin();i!=m_offsetLookup.end();++i) {
394 temp << "Tag(" << i->first << ")" << endl;
395 tempView.setOffset(i->second);
396 temp << tempView.toString() << endl;
397 }
398 return temp.str();
399 }
400
401 DataTypes::ValueType::size_type
402 DataTagged::getPointOffset(int sampleNo,
403 int dataPointNo) const
404 {
405 int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
406 DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));
407 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
408 if (pos!=m_offsetLookup.end()) {
409 offset=pos->second;
410 }
411 return offset;
412 }
413
414 DataArrayView
415 DataTagged::getDataPointByTag(int tag) const
416 {
417 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
418 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
419 if (pos!=m_offsetLookup.end()) {
420 offset=pos->second;
421 }
422 DataArrayView temp(getPointDataView());
423 temp.setOffset(offset);
424 return temp;
425 }
426
427
428 DataTypes::ValueType::const_reference
429 DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i) const
430 {
431 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
432 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
433 if (pos!=m_offsetLookup.end()) {
434 offset=pos->second;
435 }
436 DataArrayView temp(getPointDataView());
437 temp.setOffset(offset);
438 return temp.getData()[offset+i];
439 }
440
441
442 DataTypes::ValueType::reference
443 DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i)
444 {
445 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
446 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
447 if (pos!=m_offsetLookup.end()) {
448 offset=pos->second;
449 }
450 DataArrayView temp(getPointDataView());
451 temp.setOffset(offset);
452 return temp.getData()[offset+i];
453 }
454
455
456
457
458
459
460 DataArrayView
461 DataTagged::getDataPoint(int sampleNo,
462 int dataPointNo)
463 {
464 EsysAssert(validSampleNo(sampleNo),"(getDataPoint) Invalid sampleNo: " << sampleNo);
465 int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
466 return getDataPointByTag(tagKey);
467 }
468
469 int
470 DataTagged::archiveData(ofstream& archiveFile,
471 const DataTypes::ValueType::size_type noValues) const
472 {
473 return(m_data.archiveData(archiveFile, noValues));
474 }
475
476 int
477 DataTagged::extractData(ifstream& archiveFile,
478 const DataTypes::ValueType::size_type noValues)
479 {
480 return(m_data.extractData(archiveFile, noValues));
481 }
482 void
483 DataTagged::symmetric(DataAbstract* ev)
484 {
485 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
486 if (temp_ev==0) {
487 throw DataException("Error - DataTagged::symmetric casting to DataTagged failed (probably a programming error).");
488 }
489 const DataTagged::DataMapType& thisLookup=getTagLookup();
490 DataTagged::DataMapType::const_iterator i;
491 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
492 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
493 temp_ev->addTag(i->first);
494 DataArrayView thisView=getDataPointByTag(i->first);
495 DataArrayView evView=temp_ev->getDataPointByTag(i->first);
496 DataArrayView::symmetric(thisView,0,evView,0);
497 }
498 DataArrayView::symmetric(getDefaultValue(),0,temp_ev->getDefaultValue(),0);
499 }
500 void
501 DataTagged::nonsymmetric(DataAbstract* ev)
502 {
503 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
504 if (temp_ev==0) {
505 throw DataException("Error - DataTagged::nonsymmetric casting to DataTagged failed (probably a programming error).");
506 }
507 const DataTagged::DataMapType& thisLookup=getTagLookup();
508 DataTagged::DataMapType::const_iterator i;
509 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
510 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
511 temp_ev->addTag(i->first);
512 DataArrayView thisView=getDataPointByTag(i->first);
513 DataArrayView evView=temp_ev->getDataPointByTag(i->first);
514 DataArrayView::nonsymmetric(thisView,0,evView,0);
515 }
516 DataArrayView::nonsymmetric(getDefaultValue(),0,temp_ev->getDefaultValue(),0);
517 }
518 void
519 DataTagged::trace(DataAbstract* ev, int axis_offset)
520 {
521 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
522 if (temp_ev==0) {
523 throw DataException("Error - DataTagged::trace casting to DataTagged failed (probably a programming error).");
524 }
525 const DataTagged::DataMapType& thisLookup=getTagLookup();
526 DataTagged::DataMapType::const_iterator i;
527 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
528 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
529 temp_ev->addTag(i->first);
530 DataArrayView thisView=getDataPointByTag(i->first);
531 DataArrayView evView=temp_ev->getDataPointByTag(i->first);
532 DataArrayView::trace(thisView,0,evView,0, axis_offset);
533 }
534 DataArrayView::trace(getDefaultValue(),0,temp_ev->getDefaultValue(),0,axis_offset);
535 }
536
537 void
538 DataTagged::transpose(DataAbstract* ev, int axis_offset)
539 {
540 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
541 if (temp_ev==0) {
542 throw DataException("Error - DataTagged::transpose casting to DataTagged failed (probably a programming error).");
543 }
544 const DataTagged::DataMapType& thisLookup=getTagLookup();
545 DataTagged::DataMapType::const_iterator i;
546 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
547 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
548 temp_ev->addTag(i->first);
549 DataArrayView thisView=getDataPointByTag(i->first);
550 DataArrayView evView=temp_ev->getDataPointByTag(i->first);
551 DataArrayView::transpose(thisView,0,evView,0, axis_offset);
552 }
553 DataArrayView::transpose(getDefaultValue(),0,temp_ev->getDefaultValue(),0,axis_offset);
554 }
555
556 void
557 DataTagged::swapaxes(DataAbstract* ev, int axis0, int axis1)
558 {
559 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
560 if (temp_ev==0) {
561 throw DataException("Error - DataTagged::swapaxes casting to DataTagged failed (probably a programming error).");
562 }
563 const DataTagged::DataMapType& thisLookup=getTagLookup();
564 DataTagged::DataMapType::const_iterator i;
565 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
566 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
567 temp_ev->addTag(i->first);
568 DataArrayView thisView=getDataPointByTag(i->first);
569 DataArrayView evView=temp_ev->getDataPointByTag(i->first);
570 DataArrayView::swapaxes(thisView,0,evView,0,axis0,axis1);
571 }
572 DataArrayView::swapaxes(getDefaultValue(),0,temp_ev->getDefaultValue(),0,axis0,axis1);
573 }
574
575 void
576 DataTagged::eigenvalues(DataAbstract* ev)
577 {
578 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
579 if (temp_ev==0) {
580 throw DataException("Error - DataTagged::eigenvalues casting to DataTagged failed (propably a programming error).");
581 }
582 const DataTagged::DataMapType& thisLookup=getTagLookup();
583 DataTagged::DataMapType::const_iterator i;
584 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
585 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
586 temp_ev->addTag(i->first);
587 DataArrayView thisView=getDataPointByTag(i->first);
588 DataArrayView evView=temp_ev->getDataPointByTag(i->first);
589 DataArrayView::eigenvalues(thisView,0,evView,0);
590 }
591 DataArrayView::eigenvalues(getDefaultValue(),0,temp_ev->getDefaultValue(),0);
592 }
593 void
594 DataTagged::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
595 {
596 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
597 if (temp_ev==0) {
598 throw DataException("Error - DataTagged::eigenvalues_and_eigenvectors casting to DataTagged failed (propably a programming error).");
599 }
600 DataTagged* temp_V=dynamic_cast<DataTagged*>(V);
601 if (temp_V==0) {
602 throw DataException("Error - DataTagged::eigenvalues_and_eigenvectors casting to DataTagged failed (propably a programming error).");
603 }
604 const DataTagged::DataMapType& thisLookup=getTagLookup();
605 DataTagged::DataMapType::const_iterator i;
606 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
607 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
608 temp_ev->addTag(i->first);
609 temp_V->addTag(i->first);
610 DataArrayView thisView=getDataPointByTag(i->first);
611 DataArrayView evView=temp_ev->getDataPointByTag(i->first);
612 DataArrayView VView=temp_V->getDataPointByTag(i->first);
613 DataArrayView::eigenvalues_and_eigenvectors(thisView,0,evView,0,VView,0,tol);
614 }
615 DataArrayView::eigenvalues_and_eigenvectors(getDefaultValue(),0,
616 temp_ev->getDefaultValue(),0,
617 temp_V->getDefaultValue(),0,
618 tol);
619
620
621 }
622
623 void
624 DataTagged::setToZero(){
625 DataTypes::ValueType::size_type n=m_data.size();
626 for (int i=0; i<n ;++i) m_data[i]=0.;
627 }
628
629 void
630 DataTagged::dump(const std::string fileName) const
631 {
632 #ifdef PASO_MPI
633 throw DataException("Error - DataTagged:: dump is not implemented for MPI yet.");
634 #endif
635 #ifdef USE_NETCDF
636 const int ldims=DataTypes::maxRank+1;
637 const NcDim* ncdims[ldims];
638 NcVar *var, *tags_var;
639 int rank = getPointDataView().getRank();
640 int type= getFunctionSpace().getTypeCode();
641 int ndims =0;
642 long dims[ldims];
643 const double* d_ptr=&(m_data[0]);
644 DataTypes::ShapeType shape = getPointDataView().getShape();
645
646 // netCDF error handler
647 NcError err(NcError::verbose_nonfatal);
648 // Create the file.
649 NcFile dataFile(fileName.c_str(), NcFile::Replace);
650 // check if writing was successful
651 if (!dataFile.is_valid())
652 throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");
653 if (!dataFile.add_att("type_id",1) )
654 throw DataException("Error - DataTagged:: appending data type to netCDF file failed.");
655 if (!dataFile.add_att("rank",rank) )
656 throw DataException("Error - DataTagged:: appending rank attribute to netCDF file failed.");
657 if (!dataFile.add_att("function_space_type",type))
658 throw DataException("Error - DataTagged:: appending function space attribute to netCDF file failed.");
659 ndims=rank+1;
660 if ( rank >0 ) {
661 dims[0]=shape[0];
662 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
663 throw DataException("Error - DataTagged:: appending ncdimsion 0 to netCDF file failed.");
664 }
665 if ( rank >1 ) {
666 dims[1]=shape[1];
667 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
668 throw DataException("Error - DataTagged:: appending ncdimsion 1 to netCDF file failed.");
669 }
670 if ( rank >2 ) {
671 dims[2]=shape[2];
672 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
673 throw DataException("Error - DataTagged:: appending ncdimsion 2 to netCDF file failed.");
674 }
675 if ( rank >3 ) {
676 dims[3]=shape[3];
677 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
678 throw DataException("Error - DataTagged:: appending ncdimsion 3 to netCDF file failed.");
679 }
680 const DataTagged::DataMapType& thisLookup=getTagLookup();
681 DataTagged::DataMapType::const_iterator i;
682 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
683 int ntags=1;
684 for (i=thisLookup.begin();i!=thisLookupEnd;i++) ntags++;
685 int* tags =(int*) esysUtils::malloc(ntags*sizeof(int));
686 int c=1;
687 tags[0]=-1;
688 for (i=thisLookup.begin();i!=thisLookupEnd;i++) tags[c++]=i->first;
689 dims[rank]=ntags;
690 if (! (ncdims[rank] = dataFile.add_dim("num_tags", dims[rank])) )
691 {
692 esysUtils::free(tags);
693 throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");
694 }
695 if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )
696 {
697 esysUtils::free(tags);
698 throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");
699 }
700 if (! (tags_var->put(tags,dims[rank])) )
701 {
702 esysUtils::free(tags);
703 throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");
704 }
705 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
706 {
707 esysUtils::free(tags);
708 throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");
709 }
710 if (! (var->put(d_ptr,dims)) )
711 {
712 esysUtils::free(tags);
713 throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");
714 }
715 #else
716 throw DataException("Error - DataTagged:: dump is not configured with netCDF. Please contact your installation manager.");
717 #endif
718 }
719 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26