/[escript]/trunk/weipa/src/DataVar.cpp
ViewVC logotype

Contents of /trunk/weipa/src/DataVar.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (12 months, 1 week ago) by jfenwick
File size: 23346 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2018 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include <weipa/DataVar.h>
18 #include <weipa/DomainChunk.h>
19 #include <weipa/ElementData.h>
20 #include <weipa/NodeData.h>
21 #ifndef VISIT_PLUGIN
22 #include <escript/Data.h>
23 #endif
24
25 #ifdef ESYS_HAVE_NETCDF
26 #ifdef NETCDF4
27 #include <ncVar.h>
28 #include <ncDim.h>
29 #include <escript/NCHelper.h>
30 #else
31 #include <netcdfcpp.h>
32 #endif
33 #endif
34
35 #ifdef ESYS_HAVE_SILO
36 #include <silo.h>
37 #endif
38
39 #include <numeric> // for accumulate
40 #include <iostream> // for cerr
41 #include <sstream>
42 #include <stdio.h>
43
44 using namespace std;
45
46 #ifdef NETCDF4
47 using namespace netCDF;
48 #endif
49
50 namespace weipa {
51
52 //
53 // Constructor
54 //
55 DataVar::DataVar(const string& name) :
56 initialized(false), varName(name),
57 numSamples(0), rank(0), ptsPerSample(0)
58 {
59 }
60
61 //
62 // Copy constructor
63 //
64 DataVar::DataVar(const DataVar& d) :
65 initialized(d.initialized), domain(d.domain),
66 varName(d.varName), numSamples(d.numSamples),
67 rank(d.rank), ptsPerSample(d.ptsPerSample), funcSpace(d.funcSpace),
68 centering(d.centering), shape(d.shape), sampleID(d.sampleID)
69 {
70 if (numSamples > 0) {
71 CoordArray::const_iterator it;
72 for (it = d.dataArray.begin(); it != d.dataArray.end(); it++) {
73 float* c = new float[numSamples];
74 copy(*it, (*it)+numSamples, c);
75 dataArray.push_back(c);
76 }
77 }
78 }
79
80 //
81 // Destructor
82 //
83 DataVar::~DataVar()
84 {
85 cleanup();
86 }
87
88 //
89 //
90 //
91 void DataVar::cleanup()
92 {
93 CoordArray::iterator it;
94 for (it = dataArray.begin(); it != dataArray.end(); it++)
95 delete[] *it;
96 dataArray.clear();
97 shape.clear();
98 sampleID.clear();
99 numSamples = 0;
100 initialized = false;
101 }
102
103 //
104 //
105 //
106 bool DataVar::initFromEscript(escript::Data& escriptData, const_DomainChunk_ptr dom)
107 {
108 #ifndef VISIT_PLUGIN
109 cleanup();
110
111 if (!escriptData.isConstant() && !escriptData.actsExpanded()) {
112 cerr << "WARNING: Weipa only supports constant & expanded data, "
113 << "not initializing " << varName << endl;
114 return false;
115 }
116
117 domain = dom;
118 rank = escriptData.getDataPointRank();
119 ptsPerSample = escriptData.getNumDataPointsPerSample();
120 shape = escriptData.getDataPointShape();
121 funcSpace = escriptData.getFunctionSpace().getTypeCode();
122 numSamples = escriptData.getNumSamples();
123 centering = domain->getCenteringForFunctionSpace(funcSpace);
124
125 #ifdef _DEBUG
126 cout << varName << ":\t" << numSamples << " samples, "
127 << ptsPerSample << " pts/s, rank: " << rank << endl;
128 #endif
129
130 NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
131 if (nodes == NULL)
132 return false;
133
134 meshName = nodes->getName();
135 siloMeshName = nodes->getFullSiloName();
136 initialized = true;
137
138 // no samples? Nothing more to do.
139 if (numSamples == 0)
140 return true;
141
142 const escript::DataTypes::dim_t* iPtr = escriptData.getFunctionSpace().borrowSampleReferenceIDs();
143 sampleID.insert(sampleID.end(), numSamples, 0);
144 copy(iPtr, iPtr+numSamples, sampleID.begin());
145
146 size_t dimSize = 1;
147 if (rank > 0)
148 dimSize *= shape[0];
149 if (rank > 1)
150 dimSize *= shape[1];
151 if (rank > 2) {
152 cerr << "WARNING: Rank " << rank << " data is not supported!\n";
153 initialized = false;
154 }
155
156 // special case: shape=(1,) or shape=(1,1) -> convert to scalar
157 if (dimSize==1 && rank>0) {
158 rank=0;
159 shape.clear();
160 }
161
162 if (initialized) {
163 size_t dataSize = dimSize * ptsPerSample;
164 float* tempData = new float[dataSize*numSamples];
165 float* destPtr = tempData;
166 if (escriptData.isConstant()) {
167 const escript::DataTypes::real_t* values =
168 escriptData.getDataRO();
169 for (int pointNo=0; pointNo<numSamples*ptsPerSample; pointNo++) {
170 copy(values, values+dimSize, destPtr);
171 destPtr += dimSize;
172 }
173 } else {
174 for (int sampleNo=0; sampleNo<numSamples; sampleNo++) {
175 const escript::DataTypes::real_t* values =
176 escriptData.getSampleDataRO(sampleNo);
177 copy(values, values+dataSize, destPtr);
178 destPtr += dataSize;
179 }
180 }
181
182 const float* srcPtr = tempData;
183 for (size_t i=0; i < dimSize; i++, srcPtr++) {
184 float* c = averageData(srcPtr, dimSize);
185 dataArray.push_back(c);
186 }
187 delete[] tempData;
188
189 initialized = reorderSamples();
190 }
191
192 return initialized;
193
194 #else // VISIT_PLUGIN
195 return false;
196 #endif
197 }
198
199 //
200 // Initialise with domain data
201 //
202 bool DataVar::initFromMeshData(const_DomainChunk_ptr dom, const IntVec& data,
203 int fsCode, Centering c, NodeData_ptr nodes, const IntVec& id)
204 {
205 cleanup();
206
207 domain = dom;
208 rank = 0;
209 ptsPerSample = 1;
210 centering = c;
211 sampleID = id;
212 meshName = nodes->getName();
213 siloMeshName = nodes->getFullSiloName();
214 numSamples = data.size();
215
216 if (numSamples > 0) {
217 float* c = new float[numSamples];
218 dataArray.push_back(c);
219 IntVec::const_iterator it;
220 for (it=data.begin(); it != data.end(); it++)
221 *c++ = static_cast<float>(*it);
222 }
223 initialized = true;
224
225 return initialized;
226 }
227
228 //
229 // Reads variable data from dump file
230 //
231 #ifdef NETCDF4
232 bool DataVar::initFromFile(const string& filename, const_DomainChunk_ptr dom)
233 {
234 cleanup();
235
236 #if ESYS_HAVE_NETCDF
237
238 NcFile input;
239 if (!escript::openNcFile(input, filename))
240 {
241 cerr << "Could not open input file " << filename << "." << endl;
242 return false;
243 }
244
245 NcDim dim;
246 NcGroupAtt att;
247
248 att = input.getAtt("type_id");
249 int typeID;
250 att.getValues(&typeID);
251 if (typeID != 2) {
252 cerr << "WARNING: Only expanded data supported!" << endl;
253 return false;
254 }
255
256 att = input.getAtt("rank");
257 att.getValues(&rank);
258
259 dim = input.getDim("num_data_points_per_sample");
260 ptsPerSample = dim.getSize();
261
262 att = input.getAtt("function_space_type");
263 att.getValues(&funcSpace);
264
265 centering = dom->getCenteringForFunctionSpace(funcSpace);
266
267 dim = input.getDim("num_samples");
268 numSamples = dim.getSize();
269
270 #ifdef _DEBUG
271 cout << varName << ":\t" << numSamples << " samples, "
272 << ptsPerSample << " pts/s, rank: " << rank << endl;
273 #endif
274
275 domain = dom;
276 NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
277 if (nodes == NULL) {
278 return false;
279 }
280
281 meshName = nodes->getName();
282 siloMeshName = nodes->getFullSiloName();
283 initialized = true;
284
285 size_t dimSize = 1;
286 vector<long> counts;
287
288 if (rank > 0) {
289 dim = input.getDim("d0");
290 int d = dim.getSize();
291 shape.push_back(d);
292 counts.push_back(d);
293 dimSize *= d;
294 }
295 if (rank > 1) {
296 dim = input.getDim("d1");
297 int d = dim.getSize();
298 shape.push_back(d);
299 counts.push_back(d);
300 dimSize *= d;
301 }
302 if (rank > 2) {
303 cerr << "WARNING: Rank " << rank << " data is not supported!\n";
304 initialized = false;
305 }
306
307 if (initialized && numSamples > 0) {
308 sampleID.insert(sampleID.end(), numSamples, 0);
309 NcVar var = input.getVar("id");
310 var.getVar(&sampleID[0]); // numSamples
311
312 size_t dataSize = dimSize*numSamples*ptsPerSample;
313 counts.push_back(ptsPerSample);
314 counts.push_back(numSamples);
315 float* tempData = new float[dataSize];
316 var = input.getVar("data");
317 var.getVar(tempData); // &counts[0]
318
319 const float* srcPtr = tempData;
320 for (size_t i=0; i < dimSize; i++, srcPtr++) {
321 float* c = averageData(srcPtr, dimSize);
322 dataArray.push_back(c);
323 }
324 delete[] tempData;
325
326 initialized = reorderSamples();
327 }
328 #endif // ESYS_HAVE_NETCDF
329
330 return initialized;
331 }
332
333 #else
334
335 bool DataVar::initFromFile(const string& filename, const_DomainChunk_ptr dom)
336 {
337 cleanup();
338
339 #if ESYS_HAVE_NETCDF
340 NcError ncerr(NcError::silent_nonfatal);
341 NcFile* input = new NcFile(filename.c_str());
342 if (!input->is_valid()) {
343 cerr << "Could not open input file " << filename << "." << endl;
344 delete input;
345 return false;
346 }
347
348 NcDim* dim;
349 NcAtt* att;
350
351 att = input->get_att("type_id");
352 int typeID = att->as_int(0);
353 if (typeID != 2) {
354 cerr << "WARNING: Only expanded data supported!" << endl;
355 delete input;
356 return false;
357 }
358
359 att = input->get_att("rank");
360 rank = att->as_int(0);
361
362 dim = input->get_dim("num_data_points_per_sample");
363 ptsPerSample = dim->size();
364
365 att = input->get_att("function_space_type");
366 funcSpace = att->as_int(0);
367
368 centering = dom->getCenteringForFunctionSpace(funcSpace);
369
370 dim = input->get_dim("num_samples");
371 numSamples = dim->size();
372
373 #ifdef _DEBUG
374 cout << varName << ":\t" << numSamples << " samples, "
375 << ptsPerSample << " pts/s, rank: " << rank << endl;
376 #endif
377
378 domain = dom;
379 NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
380 if (nodes == NULL) {
381 delete input;
382 return false;
383 }
384
385 meshName = nodes->getName();
386 siloMeshName = nodes->getFullSiloName();
387 initialized = true;
388
389 size_t dimSize = 1;
390 vector<long> counts;
391
392 if (rank > 0) {
393 dim = input->get_dim("d0");
394 int d = dim->size();
395 shape.push_back(d);
396 counts.push_back(d);
397 dimSize *= d;
398 }
399 if (rank > 1) {
400 dim = input->get_dim("d1");
401 int d = dim->size();
402 shape.push_back(d);
403 counts.push_back(d);
404 dimSize *= d;
405 }
406 if (rank > 2) {
407 cerr << "WARNING: Rank " << rank << " data is not supported!\n";
408 initialized = false;
409 }
410
411 if (initialized && numSamples > 0) {
412 sampleID.insert(sampleID.end(), numSamples, 0);
413 NcVar* var = input->get_var("id");
414 var->get(&sampleID[0], numSamples);
415
416 size_t dataSize = dimSize*numSamples*ptsPerSample;
417 counts.push_back(ptsPerSample);
418 counts.push_back(numSamples);
419 float* tempData = new float[dataSize];
420 var = input->get_var("data");
421 var->get(tempData, &counts[0]);
422
423 const float* srcPtr = tempData;
424 for (size_t i=0; i < dimSize; i++, srcPtr++) {
425 float* c = averageData(srcPtr, dimSize);
426 dataArray.push_back(c);
427 }
428 delete[] tempData;
429
430 initialized = reorderSamples();
431 }
432
433 delete input;
434 #endif // ESYS_HAVE_NETCDF
435
436 return initialized;
437 }
438 #endif
439
440 //
441 // Returns true if the data values are nodal, false if they are zonal.
442 //
443 bool DataVar::isNodeCentered() const
444 {
445 return (centering == NODE_CENTERED);
446 }
447
448 //
449 // Returns a subset of the src array according to stride parameter.
450 // If samples consist of multiple values they are averaged beforehand.
451 // Used to separate (x0,y0,z0,x1,y1,z1,...) into (x0,x1,...), (y0,y1,...) and
452 // (z0,z1,...)
453 //
454 float* DataVar::averageData(const float* src, size_t stride)
455 {
456 float* res;
457
458 if (ptsPerSample == 1) {
459 res = new float[numSamples];
460 float* dest = res;
461 for (int i=0; i<numSamples; i++, src+=stride)
462 *dest++ = *src;
463 } else {
464 ElementData_ptr cells = domain->getElementsForFunctionSpace(funcSpace);
465 int cellFactor = cells->getElementFactor();
466 res = new float[cellFactor * numSamples];
467 float* dest = res;
468 QuadMaskInfo qmi = cells->getQuadMask(funcSpace);
469 if (!qmi.mask.empty()) {
470 const float* tmpSrc = src;
471 for (int i=0; i<numSamples; i++, tmpSrc+=stride*ptsPerSample) {
472 for (int l=0; l<cellFactor; l++) {
473 double tmpVal = 0.0;
474 for (int j=0; j<ptsPerSample; j++) {
475 if (qmi.mask[l][j] != 0) {
476 tmpVal += *(tmpSrc+stride*j);
477 }
478 }
479 *dest++ = (float)(tmpVal / qmi.factor[l]);
480 }
481 }
482 } else {
483 for (int i=0; i<numSamples; i++) {
484 double tmpVal = 0.0;
485 for (int j=0; j<ptsPerSample; j++, src+=stride) {
486 tmpVal += *src;
487 }
488 tmpVal /= ptsPerSample;
489 for (int l=0; l<cellFactor; l++) {
490 *dest++ = static_cast<float>(tmpVal);
491 }
492 }
493 }
494 }
495 return res;
496 }
497
498 //
499 // Filters and reorders the raw sample values according to the node/element
500 // IDs. This is used to have data arrays ordered according to the underlying
501 // mesh (i.e. DataID[i]==MeshNodeID[i])
502 //
503 bool DataVar::reorderSamples()
504 {
505 if (numSamples == 0)
506 return true;
507
508 const IntVec* requiredIDs = NULL;
509 int requiredNumSamples = 0;
510 int cellFactor = 1;
511
512 if (centering == NODE_CENTERED) {
513 NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
514 requiredIDs = &nodes->getNodeIDs();
515 requiredNumSamples = nodes->getNumNodes();
516 } else {
517 ElementData_ptr cells = domain->getElementsForFunctionSpace(funcSpace);
518 if (cells == NULL)
519 return false;
520
521 requiredIDs = &cells->getIDs();
522 requiredNumSamples = cells->getNumElements();
523 cellFactor = cells->getElementFactor();
524 if (cellFactor > 1) {
525 numSamples *= cellFactor;
526 // update sample IDs
527 IntVec newSampleID(numSamples);
528 IntVec::const_iterator idIt = sampleID.begin();
529 IntVec::iterator newIDit = newSampleID.begin();
530 for (; idIt != sampleID.end(); idIt++, newIDit+=cellFactor) {
531 fill(newIDit, newIDit+cellFactor, *idIt);
532 }
533 sampleID.swap(newSampleID);
534 }
535 }
536
537 if (requiredNumSamples > numSamples) {
538 cerr << "ERROR: " << varName << " has " << numSamples
539 << " instead of " << requiredNumSamples << " samples!" << endl;
540 return false;
541 }
542
543 IndexMap sampleID2idx = buildIndexMap();
544 numSamples = requiredNumSamples;
545
546 // now filter the data
547 for (size_t i=0; i < dataArray.size(); i++) {
548 float* c = new float[numSamples];
549 const float* src = dataArray[i];
550 IntVec::const_iterator idIt = requiredIDs->begin();
551 size_t destIdx = 0;
552 for (; idIt != requiredIDs->end(); idIt+=cellFactor, destIdx+=cellFactor) {
553 size_t srcIdx = sampleID2idx.find(*idIt)->second;
554 copy(&src[srcIdx], &src[srcIdx+cellFactor], &c[destIdx]);
555 }
556 delete[] dataArray[i];
557 dataArray[i] = c;
558 }
559
560 // sample IDs now = mesh node/element IDs
561 sampleID = *requiredIDs;
562
563 return true;
564 }
565
566 //
567 //
568 //
569 int DataVar::getNumberOfComponents() const
570 {
571 return (rank == 0 ? 1 : accumulate(shape.begin(), shape.end(), 0));
572 }
573
574 //
575 //
576 //
577 float* DataVar::getDataFlat() const
578 {
579 int totalSize = numSamples * getNumberOfComponents();
580 float* res = new float[totalSize];
581 if (rank == 0) {
582 copy(dataArray[0], dataArray[0]+numSamples, res);
583 } else if (rank == 1) {
584 float *dest = res;
585 for (size_t c=0; c<numSamples; c++) {
586 for (size_t i=0; i<shape[0]; i++) {
587 *dest++ = dataArray[i][c];
588 }
589 }
590 } else if (rank == 2) {
591 float *dest = res;
592 for (size_t c=0; c<numSamples; c++) {
593 for (int i=0; i<shape[1]; i++) {
594 for (int j=0; j<shape[0]; j++) {
595 *dest++ = dataArray[i*shape[0]+j][c];
596 }
597 }
598 }
599 }
600
601 return res;
602 }
603
604 //
605 //
606 //
607 void DataVar::sampleToStream(ostream& os, int index)
608 {
609 // index is -1 for dummy samples, i.e. if writing the full mesh but
610 // only a reduced number of samples is required
611 if (rank == 0) {
612 if (index < 0) {
613 os << 0.;
614 } else {
615 os << dataArray[0][index];
616 }
617 } else if (rank == 1) {
618 if (index < 0) {
619 os << 0. << " " << 0. << " " << 0.;
620 } else if (shape[0] < 3) {
621 os << dataArray[0][index] << " " << dataArray[1][index]
622 << " " << 0.;
623 } else {
624 os << dataArray[0][index] << " " << dataArray[1][index]
625 << " " << dataArray[2][index];
626 }
627 } else if (rank == 2) {
628 if (index < 0) {
629 os << 0. << " " << 0. << " " << 0. << " ";
630 os << 0. << " " << 0. << " " << 0. << " ";
631 os << 0. << " " << 0. << " " << 0.;
632 } else if (shape[1] < 3) {
633 os << dataArray[0][index] << " " << dataArray[1][index]
634 << " " << 0. << " ";
635 os << dataArray[2][index] << " " << dataArray[3][index]
636 << " " << 0. << " ";
637 os << 0. << " " << 0. << " " << 0.;
638 } else {
639 os << dataArray[0][index] << " " << dataArray[1][index]
640 << " " << dataArray[2][index] << " ";
641 os << dataArray[3][index] << " " << dataArray[4][index]
642 << " " << dataArray[5][index] << " ";
643 os << dataArray[6][index] << " " << dataArray[7][index]
644 << " " << dataArray[8][index];
645 }
646 }
647 os << endl;
648 }
649
650 //
651 //
652 //
653 void DataVar::writeToVTK(ostream& os, int ownIndex)
654 {
655 if (numSamples == 0)
656 return;
657
658 if (isNodeCentered()) {
659 // data was reordered in reorderSamples() but for VTK we write the
660 // original node mesh and thus need the original ordering.
661 // Note, that this also means we may not have samples for all nodes
662 // in which case we set idx to -1 and write a dummy sample
663 const IntVec& requiredIDs = domain->getNodes()->getNodeIDs();
664 const IntVec& nodeGNI = domain->getNodes()->getGlobalNodeIndices();
665 const IntVec& nodeDist = domain->getNodes()->getNodeDistribution();
666 int firstId = nodeDist[ownIndex];
667 int lastId = nodeDist[ownIndex+1];
668 IndexMap sampleID2idx = buildIndexMap();
669 for (size_t i=0; i<nodeGNI.size(); i++) {
670 if (firstId <= nodeGNI[i] && nodeGNI[i] < lastId) {
671 IndexMap::const_iterator it = sampleID2idx.find(requiredIDs[i]);
672 int idx = (it==sampleID2idx.end() ? -1 : (int)it->second);
673 sampleToStream(os, idx);
674 }
675 }
676 } else {
677 // cell data: ghost cells have been removed so do not write ghost
678 // samples (which are the last elements in the arrays)
679 int toWrite = domain->getElementsByName(meshName)->getNumElements();
680 for (int i=0; i<toWrite; i++) {
681 sampleToStream(os, i);
682 }
683 }
684 }
685
686 ///////////////////////////////
687 // SILO related methods follow
688 ///////////////////////////////
689
690 //
691 // If the data is tensor data then the components of the tensor are stored
692 // separately in the Silo file. This method then returns a string that
693 // contains the proper Silo expression to put the tensor together again.
694 // For non-tensor data this method returns an empty string.
695 //
696 string DataVar::getTensorDef() const
697 {
698 if (rank < 2 || !initialized)
699 return string();
700
701 /// Format string for Silo 2x2 tensor
702 const string tensor2DefFmt =
703 "{{ <%sa_00>, <%sa_01> },"
704 " { <%sa_10>, <%sa_11> }}";
705
706 /// Format string for Silo 3x3 tensor
707 const string tensor3DefFmt =
708 "{{ <%sa_00>, <%sa_01>, <%sa_02> },"
709 " { <%sa_10>, <%sa_11>, <%sa_12> },"
710 " { <%sa_20>, <%sa_21>, <%sa_22> }}";
711
712 string tensorDef;
713 string tensorDir = varName+string("_comps/");
714 if (shape[1] == 3) {
715 char* tDef = new char[tensor3DefFmt.length()+9*tensorDir.length()];
716 sprintf(tDef, tensor3DefFmt.c_str(),
717 tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
718 tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str(),
719 tensorDir.c_str(), tensorDir.c_str(), tensorDir.c_str());
720 tensorDef = tDef;
721 delete[] tDef;
722 } else {
723 char* tDef = new char[tensor2DefFmt.length()+4*tensorDir.length()];
724 sprintf(tDef, tensor2DefFmt.c_str(),
725 tensorDir.c_str(), tensorDir.c_str(),
726 tensorDir.c_str(), tensorDir.c_str(),
727 tensorDir.c_str(), tensorDir.c_str());
728 tensorDef = tDef;
729 delete[] tDef;
730 }
731 return tensorDef;
732 }
733
734 //
735 // Writes the data to given Silo file under the virtual path provided.
736 // The corresponding mesh must have been written already.
737 //
738 bool DataVar::writeToSilo(DBfile* dbfile, const string& siloPath,
739 const string& units)
740 {
741 #ifdef ESYS_HAVE_SILO
742 if (!initialized)
743 return false;
744
745 if (numSamples == 0)
746 return true;
747
748 int ret;
749
750 if (siloPath != "") {
751 ret = DBSetDir(dbfile, siloPath.c_str());
752 if (ret != 0)
753 return false;
754 }
755
756 char* siloMesh = const_cast<char*>(siloMeshName.c_str());
757 int dcenter = (centering == NODE_CENTERED ? DB_NODECENT : DB_ZONECENT);
758 DBoptlist* optList = DBMakeOptlist(2);
759 if (units.length()>0) {
760 DBAddOption(optList, DBOPT_UNITS, (void*)units.c_str());
761 }
762
763 if (rank == 0) {
764 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMesh, dataArray[0],
765 numSamples, NULL, 0, DB_FLOAT, dcenter, optList);
766 }
767 else if (rank == 1) {
768 const string comps[3] = {
769 varName+string("_x"), varName+string("_y"), varName+string("_z")
770 };
771 const char* varnames[3] = {
772 comps[0].c_str(), comps[1].c_str(), comps[2].c_str()
773 };
774
775 ret = DBPutUcdvar(dbfile, varName.c_str(), siloMesh, shape[0],
776 (char**)varnames, &dataArray[0], numSamples, NULL,
777 0, DB_FLOAT, dcenter, optList);
778 }
779 else {
780 string tensorDir = varName+string("_comps/");
781 ret = DBMkdir(dbfile, tensorDir.c_str());
782 if (ret == 0) {
783 int one = 1;
784 DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
785
786 for (int i=0; i<shape[1]; i++) {
787 for (int j=0; j<shape[0]; j++) {
788 ostringstream varname;
789 varname << tensorDir << "a_" << i << j;
790 ret = DBPutUcdvar1(dbfile, varname.str().c_str(), siloMesh,
791 dataArray[i*shape[0]+j], numSamples,
792 NULL, 0, DB_FLOAT, dcenter, optList);
793 if (ret != 0) break;
794 }
795 if (ret != 0) break;
796 }
797 } // ret==0
798 } // rank
799
800 DBFreeOptlist(optList);
801 DBSetDir(dbfile, "/");
802 return (ret == 0);
803
804 #else // !ESYS_HAVE_SILO
805 return false;
806 #endif
807 }
808
809 } // namespace weipa
810

  ViewVC Help
Powered by ViewVC 1.1.26