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 |
// ElementData.cpp |
16 |
// |
17 |
#include <escriptreader/ElementData.h> |
18 |
#include <netcdf.hh> |
19 |
#if HAVE_SILO |
20 |
#include <silo.h> |
21 |
#endif |
22 |
|
23 |
using namespace std; |
24 |
|
25 |
// The following arrays contain indices to convert unsupported element |
26 |
// types into supported ones (lines, triangles, quads, hexahedrons) |
27 |
|
28 |
static const size_t line3indices[2*2] = { |
29 |
0, 2, |
30 |
2, 1 |
31 |
}; |
32 |
static const size_t tri6indices[4*3] = { |
33 |
0, 3, 5, |
34 |
3, 1, 4, |
35 |
5, 4, 2, |
36 |
3, 4, 5 |
37 |
}; |
38 |
static const size_t rec8indices[6*3] = { |
39 |
0, 4, 7, |
40 |
4, 1, 5, |
41 |
5, 2, 6, |
42 |
6, 3, 7, |
43 |
4, 6, 7, |
44 |
4, 5, 6 |
45 |
}; |
46 |
static const size_t rec9indices[4*4] = { |
47 |
0, 4, 8, 7, |
48 |
4, 1, 5, 8, |
49 |
7, 8, 6, 3, |
50 |
8, 5, 2, 6 |
51 |
}; |
52 |
static const size_t tet10indices[8*4] = { |
53 |
6, 0, 7, 4, |
54 |
2, 6, 9, 5, |
55 |
9, 7, 3, 8, |
56 |
4, 1, 8, 5, |
57 |
6, 7, 9, 8, |
58 |
5, 6, 9, 8, |
59 |
6, 5, 4, 8, |
60 |
7, 6, 4, 8 |
61 |
}; |
62 |
static const size_t hex20indices[36*3] = { |
63 |
0, 8, 12, 8, 1, 13, 13, 5, 16, |
64 |
16, 4, 12, 8, 13, 16, 8, 16, 12, |
65 |
1, 9, 13, 9, 2, 14, 14, 6, 17, |
66 |
17, 5, 13, 9, 14, 17, 9, 17, 13, |
67 |
2, 10, 14, 10, 3, 15, 15, 7, 18, |
68 |
18, 14, 6, 10, 15, 18, 10, 18, 14, |
69 |
3, 11, 15, 11, 0, 12, 12, 4, 19, |
70 |
19, 7, 15, 11, 12, 19, 11, 19, 15, |
71 |
4, 16, 19, 16, 5, 17, 17, 6, 18, |
72 |
18, 7, 19, 16, 17, 18, 16, 18, 19, |
73 |
3, 10, 11, 10, 2, 9, 9, 1, 8, |
74 |
8, 0, 11, 10, 9, 8, 10, 8, 11 |
75 |
}; |
76 |
static const size_t hex27indices[8*8] = { |
77 |
0, 8, 20, 11, 12, 21, 26, 24, |
78 |
8, 1, 9, 20, 21, 13, 22, 26, |
79 |
11, 20, 10, 3, 24, 26, 23, 15, |
80 |
20, 9, 2, 10, 26, 22, 14, 23, |
81 |
12, 21, 26, 24, 4, 16, 25, 19, |
82 |
21, 13, 22, 26, 16, 5, 17, 25, |
83 |
24, 26, 23, 15, 19, 25, 18, 7, |
84 |
26, 22, 14, 23, 25, 17, 6, 18 |
85 |
}; |
86 |
|
87 |
namespace EscriptReader { |
88 |
|
89 |
// |
90 |
// Constructor |
91 |
// |
92 |
ElementData::ElementData(const string& elementName, const Mesh* mainMesh) |
93 |
: name(elementName), count(0), reducedCount(0), numGhostElements(0), |
94 |
numReducedGhostElements(0), fullMesh(NULL), reducedMesh(NULL), |
95 |
originalMesh(mainMesh), fullMeshIsOriginalMesh(false) |
96 |
{ |
97 |
} |
98 |
|
99 |
// |
100 |
// Copy constructor |
101 |
// |
102 |
ElementData::ElementData(const ElementData& e) |
103 |
{ |
104 |
name = e.name; |
105 |
count = e.count; |
106 |
reducedCount = e.reducedCount; |
107 |
numGhostElements = e.numGhostElements; |
108 |
numReducedGhostElements = e.numReducedGhostElements; |
109 |
numDims = e.numDims; |
110 |
type = e.type; |
111 |
reducedType = e.reducedType; |
112 |
nodesPerElement = e.nodesPerElement; |
113 |
reducedNodesPerElement = e.reducedNodesPerElement; |
114 |
|
115 |
originalMesh = e.originalMesh; |
116 |
fullMeshIsOriginalMesh = e.fullMeshIsOriginalMesh; |
117 |
if (fullMeshIsOriginalMesh) |
118 |
fullMesh = const_cast<Mesh*>(originalMesh); |
119 |
else if (e.fullMesh) |
120 |
fullMesh = new Mesh(*e.fullMesh); |
121 |
else |
122 |
fullMesh = NULL; |
123 |
|
124 |
if (e.reducedMesh) |
125 |
reducedMesh = new Mesh(*e.reducedMesh); |
126 |
else |
127 |
reducedMesh = NULL; |
128 |
|
129 |
nodes = e.nodes; |
130 |
reducedNodes = e.reducedNodes; |
131 |
ID = e.ID; |
132 |
color = e.color; |
133 |
tag = e.tag; |
134 |
owner = e.owner; |
135 |
reducedOwner = e.reducedOwner; |
136 |
indexArray = e.indexArray; |
137 |
reducedIndexArray = e.reducedIndexArray; |
138 |
ID2idx = e.ID2idx; |
139 |
} |
140 |
|
141 |
// |
142 |
// Destructor |
143 |
// |
144 |
ElementData::~ElementData() |
145 |
{ |
146 |
if (reducedMesh) |
147 |
delete reducedMesh; |
148 |
if (fullMesh && !fullMeshIsOriginalMesh) |
149 |
delete fullMesh; |
150 |
} |
151 |
|
152 |
StringVec ElementData::getMeshNames() const |
153 |
{ |
154 |
StringVec res; |
155 |
if (fullMesh && !fullMeshIsOriginalMesh) |
156 |
res.push_back(fullMesh->getName()); |
157 |
if (reducedMesh) |
158 |
res.push_back(reducedMesh->getName()); |
159 |
return res; |
160 |
} |
161 |
|
162 |
StringVec ElementData::getVarNames() const |
163 |
{ |
164 |
StringVec res; |
165 |
if (count > 0) { |
166 |
res.push_back(name + string("_Color")); |
167 |
res.push_back(name + string("_Id")); |
168 |
res.push_back(name + string("_Owner")); |
169 |
res.push_back(name + string("_Tag")); |
170 |
} |
171 |
return res; |
172 |
} |
173 |
|
174 |
const IntVec& ElementData::getVarDataByName(const string varName) const |
175 |
{ |
176 |
if (varName == name+string("_Color")) |
177 |
return color; |
178 |
else if (varName == name+string("_Id")) |
179 |
return ID; |
180 |
else if (varName == name+string("_Owner")) { |
181 |
if (reducedCount > 0) |
182 |
return reducedOwner; |
183 |
else |
184 |
return owner; |
185 |
} else if (varName == name+string("_Tag")) |
186 |
return tag; |
187 |
else |
188 |
throw "Invalid variable name"; |
189 |
} |
190 |
|
191 |
// |
192 |
// Reads element data from given NetCDF file |
193 |
// |
194 |
bool ElementData::readFromNc(NcFile* ncfile) |
195 |
{ |
196 |
NcAtt* att; |
197 |
NcVar* var; |
198 |
|
199 |
string num_str("num_"); |
200 |
num_str += name; |
201 |
|
202 |
att = ncfile->get_att(num_str.c_str()); |
203 |
count = att->as_int(0); |
204 |
|
205 |
// Only attempt to read data if there are any elements. |
206 |
// Having no elements is not an error. |
207 |
if (count == 0) |
208 |
return true; |
209 |
|
210 |
att = ncfile->get_att((num_str + string("_numNodes")).c_str()); |
211 |
nodesPerElement = att->as_int(0); |
212 |
|
213 |
nodes.insert(nodes.end(), count*nodesPerElement, 0); |
214 |
var = ncfile->get_var((name + string("_Nodes")).c_str()); |
215 |
var->get(&nodes[0], count, nodesPerElement); |
216 |
|
217 |
color.insert(color.end(), count, 0); |
218 |
var = ncfile->get_var((name + string("_Color")).c_str()); |
219 |
var->get(&color[0], count); |
220 |
|
221 |
ID.insert(ID.end(), count, 0); |
222 |
var = ncfile->get_var((name + string("_Id")).c_str()); |
223 |
var->get(&ID[0], count); |
224 |
|
225 |
owner.insert(owner.end(), count, 0); |
226 |
var = ncfile->get_var((name + string("_Owner")).c_str()); |
227 |
var->get(&owner[0], count); |
228 |
|
229 |
tag.insert(tag.end(), count, 0); |
230 |
var = ncfile->get_var((name + string("_Tag")).c_str()); |
231 |
var->get(&tag[0], count); |
232 |
|
233 |
att = ncfile->get_att((name + string("_TypeId")).c_str()); |
234 |
FinleyElementInfo f = getFinleyTypeInfo((ElementTypeId)att->as_int(0)); |
235 |
type = f.elementType; |
236 |
reducedType = f.reducedElementType; |
237 |
|
238 |
// build elementID->index map |
239 |
buildIndexMap(); |
240 |
|
241 |
if (f.elementFactor > 1 || f.reducedElementSize != nodesPerElement) |
242 |
buildReducedElements(f); |
243 |
|
244 |
buildMeshes(); |
245 |
|
246 |
return true; |
247 |
} |
248 |
|
249 |
// |
250 |
// |
251 |
// |
252 |
void ElementData::reorderArray(IntVec& v, const IntVec& idx, |
253 |
int elementsPerIndex) |
254 |
{ |
255 |
IntVec newArray(v.size()); |
256 |
IntVec::iterator arrIt = newArray.begin(); |
257 |
IntVec::const_iterator idxIt; |
258 |
if (elementsPerIndex == 1) { |
259 |
for (idxIt=idx.begin(); idxIt!=idx.end(); idxIt++) { |
260 |
*arrIt++ = v[*idxIt]; |
261 |
} |
262 |
} else { |
263 |
for (idxIt=idx.begin(); idxIt!=idx.end(); idxIt++) { |
264 |
int i = *idxIt; |
265 |
copy(&v[i*elementsPerIndex], &v[(i+1)*elementsPerIndex], arrIt); |
266 |
arrIt += elementsPerIndex; |
267 |
} |
268 |
} |
269 |
v.swap(newArray); |
270 |
} |
271 |
|
272 |
// |
273 |
// |
274 |
// |
275 |
void ElementData::buildReducedElements(const FinleyElementInfo& f) |
276 |
{ |
277 |
reducedNodes.clear(); |
278 |
reducedNodes.insert(reducedNodes.end(), |
279 |
f.reducedElementSize*count, 0); |
280 |
IntVec::iterator reducedIt = reducedNodes.begin(); |
281 |
IntVec::const_iterator origIt; |
282 |
for (origIt=nodes.begin(); origIt!=nodes.end(); |
283 |
origIt+=nodesPerElement) |
284 |
{ |
285 |
std::copy(origIt, origIt+f.reducedElementSize, reducedIt); |
286 |
reducedIt += f.reducedElementSize; |
287 |
} |
288 |
|
289 |
// Remove comment to save disk space - we don't really need the full |
290 |
// elements except for the main mesh |
291 |
if (f.elementFactor > 1 /*&& name == "Elements"*/) { |
292 |
// replace each element by multiple smaller ones |
293 |
IntVec fullNodes(f.elementSize*f.elementFactor*count); |
294 |
IntVec::iterator cellIt = fullNodes.begin(); |
295 |
|
296 |
// copy owner data |
297 |
owner.swap(reducedOwner); |
298 |
owner.clear(); |
299 |
for (int i=0; i < count; i++) { |
300 |
owner.insert(owner.end(), f.elementFactor, reducedOwner[i]); |
301 |
for (int j=0; j < f.elementFactor*f.elementSize; j++) |
302 |
*cellIt++ = nodes[ |
303 |
nodesPerElement*i+f.multiCellIndices[j]]; |
304 |
} |
305 |
nodes.swap(fullNodes); |
306 |
reducedCount = count; |
307 |
reducedNodesPerElement = f.reducedElementSize; |
308 |
nodesPerElement = f.elementSize; |
309 |
count *= f.elementFactor; |
310 |
} else { |
311 |
// we only keep the reduced elements but treat them as regular |
312 |
// ones, so replace the data accordingly |
313 |
nodes.swap(reducedNodes); |
314 |
reducedNodes.clear(); |
315 |
nodesPerElement = f.reducedElementSize; |
316 |
type = f.reducedElementType; |
317 |
} |
318 |
} |
319 |
|
320 |
// |
321 |
// |
322 |
// |
323 |
void ElementData::prepareGhostIndices(int ownIndex) |
324 |
{ |
325 |
indexArray.clear(); |
326 |
reducedIndexArray.clear(); |
327 |
numGhostElements = 0; |
328 |
numReducedGhostElements = 0; |
329 |
|
330 |
// move indices of "ghost zones" to the end to be able to reorder |
331 |
// data accordingly |
332 |
for (int i=0; i<count; i++) |
333 |
if (owner[i] == ownIndex) |
334 |
indexArray.push_back(i); |
335 |
|
336 |
for (int i=0; i<count; i++) |
337 |
if (owner[i] != ownIndex) { |
338 |
numGhostElements++; |
339 |
indexArray.push_back(i); |
340 |
} |
341 |
|
342 |
for (int i=0; i<reducedCount; i++) |
343 |
if (reducedOwner[i] == ownIndex) |
344 |
reducedIndexArray.push_back(i); |
345 |
|
346 |
for (int i=0; i<reducedCount; i++) |
347 |
if (reducedOwner[i] != ownIndex) { |
348 |
numReducedGhostElements++; |
349 |
reducedIndexArray.push_back(i); |
350 |
} |
351 |
} |
352 |
|
353 |
// |
354 |
// |
355 |
// |
356 |
void ElementData::handleGhostZones(int ownIndex) |
357 |
{ |
358 |
prepareGhostIndices(ownIndex); |
359 |
|
360 |
// move "ghost data" to the end of the arrays |
361 |
if (numGhostElements > 0) { |
362 |
reorderArray(nodes, indexArray, nodesPerElement); |
363 |
reorderArray(owner, indexArray, 1); |
364 |
if (reducedCount == 0) { |
365 |
reorderArray(color, indexArray, 1); |
366 |
reorderArray(ID, indexArray, 1); |
367 |
reorderArray(tag, indexArray, 1); |
368 |
} |
369 |
} |
370 |
|
371 |
if (numReducedGhostElements > 0) { |
372 |
reorderArray(reducedNodes, reducedIndexArray, reducedNodesPerElement); |
373 |
reorderArray(reducedOwner, reducedIndexArray, 1); |
374 |
reorderArray(color, reducedIndexArray, 1); |
375 |
reorderArray(ID, reducedIndexArray, 1); |
376 |
reorderArray(tag, reducedIndexArray, 1); |
377 |
} |
378 |
} |
379 |
|
380 |
// |
381 |
// |
382 |
// |
383 |
void ElementData::removeGhostZones() |
384 |
{ |
385 |
if (numGhostElements > 0) { |
386 |
count -= numGhostElements; |
387 |
nodes.resize(count*nodesPerElement); |
388 |
owner.resize(count); |
389 |
if (reducedCount == 0) { |
390 |
color.resize(count); |
391 |
ID.resize(count); |
392 |
tag.resize(count); |
393 |
} |
394 |
numGhostElements = 0; |
395 |
} |
396 |
|
397 |
if (numReducedGhostElements > 0) { |
398 |
reducedCount -= numReducedGhostElements; |
399 |
reducedNodes.resize(reducedCount*reducedNodesPerElement); |
400 |
reducedOwner.resize(reducedCount); |
401 |
color.resize(reducedCount); |
402 |
ID.resize(reducedCount); |
403 |
tag.resize(reducedCount); |
404 |
numReducedGhostElements = 0; |
405 |
} |
406 |
} |
407 |
|
408 |
// |
409 |
// |
410 |
// |
411 |
void ElementData::buildMeshes() |
412 |
{ |
413 |
if (count == 0) |
414 |
return; |
415 |
|
416 |
// use existing original mesh for Elements but build a new mesh |
417 |
// containing only the required nodes for other types |
418 |
if (name == "Elements") { |
419 |
fullMesh = const_cast<Mesh*>(originalMesh); |
420 |
fullMeshIsOriginalMesh = true; |
421 |
} else { |
422 |
// first: build a map of required IDs while translating |
423 |
// the original node IDs at the same time |
424 |
IntVec::iterator it; |
425 |
IndexMap nodeID2idx; |
426 |
IntVec nodeIDs; |
427 |
size_t newIdx = 0; |
428 |
|
429 |
for (it = nodes.begin(); it != nodes.end(); it++) { |
430 |
IndexMap::iterator res = nodeID2idx.find(*it); |
431 |
if (res == nodeID2idx.end()) { |
432 |
nodeIDs.push_back(*it); |
433 |
nodeID2idx[*it] = newIdx; |
434 |
*it = newIdx++; |
435 |
} else { |
436 |
*it = res->second; |
437 |
} |
438 |
} |
439 |
|
440 |
// second: use map to fill coordinates for the nodes |
441 |
// newIdx contains the number of nodes |
442 |
CoordArray coords; |
443 |
const CoordArray& origCoords = originalMesh->getCoords(); |
444 |
for (int dim=0; dim < originalMesh->getNumDims(); dim++) { |
445 |
float* c = new float[newIdx]; |
446 |
coords.push_back(c); |
447 |
IndexMap::const_iterator mIt; |
448 |
for (mIt = nodeID2idx.begin(); mIt != nodeID2idx.end(); mIt++) |
449 |
c[mIt->second] = origCoords[dim][mIt->first]; |
450 |
} |
451 |
|
452 |
if (fullMesh) |
453 |
delete fullMesh; |
454 |
|
455 |
fullMesh = new Mesh(coords, originalMesh->getNumDims(), newIdx); |
456 |
fullMesh->setName(name); |
457 |
fullMesh->setIndexMap(nodeID2idx); |
458 |
fullMesh->setNodeIDs(nodeIDs); |
459 |
} |
460 |
|
461 |
#ifdef _DEBUG |
462 |
cout << fullMesh->getName() << " has " << fullMesh->getNumNodes() |
463 |
<< " nodes, " << count << " elements" << endl; |
464 |
#endif |
465 |
|
466 |
// build a reduced mesh if necessary |
467 |
if (reducedCount > 0) { |
468 |
IntVec::iterator it; |
469 |
IndexMap nodeID2idx; |
470 |
IntVec reducedNodeIDs; |
471 |
size_t newIdx = 0; |
472 |
|
473 |
for (it = reducedNodes.begin(); it != reducedNodes.end(); it++) { |
474 |
int id = originalMesh->getNodeIDs()[*it]; |
475 |
IndexMap::iterator res = nodeID2idx.find(id); |
476 |
if (res == nodeID2idx.end()) { |
477 |
reducedNodeIDs.push_back(id); |
478 |
nodeID2idx[id] = newIdx; |
479 |
*it = newIdx++; |
480 |
} else { |
481 |
*it = res->second; |
482 |
} |
483 |
} |
484 |
|
485 |
CoordArray coords; |
486 |
const CoordArray& origCoords = originalMesh->getCoords(); |
487 |
for (int dim=0; dim < originalMesh->getNumDims(); dim++) { |
488 |
float* c = new float[newIdx]; |
489 |
coords.push_back(c); |
490 |
IndexMap::const_iterator mIt; |
491 |
for (mIt = nodeID2idx.begin(); mIt != nodeID2idx.end(); mIt++) { |
492 |
IndexMap::const_iterator idx; |
493 |
idx = originalMesh->getIndexMap().find(mIt->first); |
494 |
c[mIt->second] = origCoords[dim][idx->second]; |
495 |
} |
496 |
} |
497 |
if (reducedMesh) |
498 |
delete reducedMesh; |
499 |
|
500 |
reducedMesh = new Mesh(coords, originalMesh->getNumDims(), newIdx); |
501 |
reducedMesh->setName(string("Reduced") + name); |
502 |
reducedMesh->setIndexMap(nodeID2idx); |
503 |
reducedMesh->setNodeIDs(reducedNodeIDs); |
504 |
|
505 |
#ifdef _DEBUG |
506 |
cout << reducedMesh->getName() << " has " << newIdx |
507 |
<< " nodes, " << reducedCount << " elements" << endl; |
508 |
#endif |
509 |
} |
510 |
} |
511 |
|
512 |
#if HAVE_SILO |
513 |
// |
514 |
// |
515 |
// |
516 |
inline int toSiloElementType(int type) |
517 |
{ |
518 |
switch (type) { |
519 |
case ZONETYPE_BEAM: return DB_ZONETYPE_BEAM; |
520 |
case ZONETYPE_HEX: return DB_ZONETYPE_HEX; |
521 |
case ZONETYPE_POLYGON: return DB_ZONETYPE_POLYGON; |
522 |
case ZONETYPE_QUAD: return DB_ZONETYPE_QUAD; |
523 |
case ZONETYPE_TET: return DB_ZONETYPE_TET; |
524 |
case ZONETYPE_TRIANGLE: return DB_ZONETYPE_TRIANGLE; |
525 |
} |
526 |
return 0; |
527 |
} |
528 |
#endif |
529 |
|
530 |
// |
531 |
// |
532 |
// |
533 |
bool ElementData::writeToSilo(DBfile* dbfile, const string& siloPath) |
534 |
{ |
535 |
#if HAVE_SILO |
536 |
if (count == 0) |
537 |
return true; |
538 |
|
539 |
const char* meshName; |
540 |
int numCells; |
541 |
string varName, siloMeshName; |
542 |
int ret; |
543 |
|
544 |
if (siloPath != "") { |
545 |
ret = DBSetDir(dbfile, siloPath.c_str()); |
546 |
if (ret != 0) |
547 |
return false; |
548 |
} |
549 |
|
550 |
// write out the full mesh in any case |
551 |
if (siloPath == "/") |
552 |
siloMeshName = siloPath + fullMesh->getName(); |
553 |
else |
554 |
siloMeshName = siloPath + string("/") + fullMesh->getName(); |
555 |
|
556 |
int arraylen = count * nodesPerElement; |
557 |
int eltype = toSiloElementType(type); |
558 |
varName = name + string("_zones"); |
559 |
ret = DBPutZonelist2(dbfile, varName.c_str(), count, |
560 |
originalMesh->getNumDims(), &nodes[0], arraylen, 0, 0, |
561 |
numGhostElements, &eltype, &nodesPerElement, &count, 1, NULL); |
562 |
if (ret == 0) { |
563 |
CoordArray& coordbase = const_cast<CoordArray&>(fullMesh->getCoords()); |
564 |
ret = DBPutUcdmesh(dbfile, siloMeshName.c_str(), |
565 |
originalMesh->getNumDims(), NULL, &coordbase[0], |
566 |
fullMesh->getNumNodes(), count, varName.c_str(), |
567 |
/*"facelist"*/NULL, DB_FLOAT, NULL); |
568 |
} |
569 |
|
570 |
// Point mesh is useful for debugging |
571 |
if (0) { |
572 |
CoordArray& coordbase = const_cast<CoordArray&>(fullMesh->getCoords()); |
573 |
DBPutPointmesh(dbfile, "/pointmesh", |
574 |
originalMesh->getNumDims(), &coordbase[0], |
575 |
fullMesh->getNumNodes(), DB_FLOAT, NULL); |
576 |
} |
577 |
|
578 |
if (ret != 0) |
579 |
return false; |
580 |
|
581 |
// decide whether to additionally write out the reduced mesh |
582 |
if (reducedCount > 0) { |
583 |
if (siloPath == "/") |
584 |
siloMeshName = siloPath + reducedMesh->getName(); |
585 |
else |
586 |
siloMeshName = siloPath + string("/") + reducedMesh->getName(); |
587 |
arraylen = reducedCount * reducedNodesPerElement; |
588 |
eltype = toSiloElementType(reducedType); |
589 |
varName = string("Reduced") + name + string("_zones"); |
590 |
ret = DBPutZonelist2(dbfile, varName.c_str(), reducedCount, |
591 |
originalMesh->getNumDims(), &reducedNodes[0], arraylen, 0, 0, |
592 |
numReducedGhostElements, &eltype, &reducedNodesPerElement, |
593 |
&reducedCount, 1, NULL); |
594 |
if (ret == 0) { |
595 |
CoordArray& coordbase = const_cast<CoordArray&>(reducedMesh->getCoords()); |
596 |
ret = DBPutUcdmesh(dbfile, siloMeshName.c_str(), |
597 |
originalMesh->getNumDims(), NULL, &coordbase[0], |
598 |
reducedMesh->getNumNodes(), reducedCount, varName.c_str(), |
599 |
NULL, DB_FLOAT, NULL); |
600 |
} |
601 |
if (ret != 0) |
602 |
return false; |
603 |
numCells = reducedCount; |
604 |
} else { |
605 |
numCells = count; |
606 |
} |
607 |
meshName = siloMeshName.c_str(); |
608 |
|
609 |
// finally, write out the element-centered variables on the correct mesh |
610 |
varName = name + string("_Color"); |
611 |
ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName, |
612 |
(float*)&color[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL); |
613 |
if (ret == 0) { |
614 |
varName = name + string("_Id"); |
615 |
ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName, |
616 |
(float*)&ID[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL); |
617 |
} |
618 |
if (ret == 0) { |
619 |
varName = name + string("_Owner"); |
620 |
ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName, |
621 |
(float*)&owner[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL); |
622 |
} |
623 |
if (ret == 0) { |
624 |
varName = name + string("_Tag"); |
625 |
ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName, |
626 |
(float*)&tag[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL); |
627 |
} |
628 |
|
629 |
DBSetDir(dbfile, "/"); |
630 |
return (ret == 0); |
631 |
|
632 |
#else // !HAVE_SILO |
633 |
return false; |
634 |
#endif |
635 |
} |
636 |
|
637 |
// |
638 |
// |
639 |
// |
640 |
FinleyElementInfo ElementData::getFinleyTypeInfo(ElementTypeId typeId) |
641 |
{ |
642 |
FinleyElementInfo ret; |
643 |
ret.multiCellIndices = NULL; |
644 |
ret.elementFactor = 1; |
645 |
|
646 |
switch (typeId) { |
647 |
case Point1_Contact://untested |
648 |
case Line2Face_Contact://untested |
649 |
case Line3Face_Contact://untested |
650 |
case Line2Face://untested |
651 |
case Line3Face://untested |
652 |
case Point1://untested |
653 |
cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl; |
654 |
ret.elementSize = 1; |
655 |
ret.elementType = ZONETYPE_POLYGON; |
656 |
break; |
657 |
|
658 |
case Tri3Face_Contact://untested |
659 |
case Tri3Face://untested |
660 |
cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl; |
661 |
case Line2_Contact: |
662 |
case Rec4Face_Contact: |
663 |
case Rec4Face: |
664 |
case Line2: |
665 |
ret.elementSize = ret.reducedElementSize = 2; |
666 |
ret.elementType = ret.reducedElementType = ZONETYPE_BEAM; |
667 |
break; |
668 |
|
669 |
case Line3: |
670 |
ret.multiCellIndices = line3indices; |
671 |
ret.elementFactor = 2; |
672 |
// fall through |
673 |
case Line3_Contact: |
674 |
case Tri6Face_Contact://untested |
675 |
case Rec8Face_Contact: |
676 |
case Tri6Face://untested |
677 |
case Rec8Face: |
678 |
//VTK_QUADRATIC_EDGE |
679 |
ret.elementSize = ret.reducedElementSize = 2; |
680 |
ret.elementType = ret.reducedElementType = ZONETYPE_BEAM; |
681 |
break; |
682 |
|
683 |
case Tet4Face_Contact://untested |
684 |
case Tet4Face://untested |
685 |
cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl; |
686 |
case Tri3_Contact: |
687 |
case Tri3: |
688 |
ret.elementSize = ret.reducedElementSize = 3; |
689 |
ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE; |
690 |
break; |
691 |
|
692 |
case Rec4_Contact: |
693 |
case Hex8Face_Contact: |
694 |
case Hex8Face: |
695 |
case Rec4: |
696 |
ret.elementSize = ret.reducedElementSize = 4; |
697 |
ret.elementType = ret.reducedElementType = ZONETYPE_QUAD; |
698 |
break; |
699 |
|
700 |
case Rec9: |
701 |
ret.multiCellIndices = rec9indices; |
702 |
ret.elementFactor = 4; |
703 |
// fall through |
704 |
case Rec9_Contact: |
705 |
ret.elementSize = ret.reducedElementSize = 4; |
706 |
ret.elementType = ret.reducedElementType = ZONETYPE_QUAD; |
707 |
break; |
708 |
|
709 |
case Tet4: |
710 |
ret.elementSize = ret.reducedElementSize = 4; |
711 |
ret.elementType = ret.reducedElementType = ZONETYPE_TET; |
712 |
break; |
713 |
|
714 |
case Tri6: |
715 |
ret.multiCellIndices = tri6indices; |
716 |
ret.elementFactor = 4; |
717 |
// fall through |
718 |
case Tri6_Contact: |
719 |
case Tet10Face_Contact://untested |
720 |
case Tet10Face://untested |
721 |
//VTK_QUADRATIC_TRIANGLE |
722 |
ret.elementSize = ret.reducedElementSize = 3; |
723 |
ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE; |
724 |
break; |
725 |
|
726 |
case Rec8: |
727 |
ret.multiCellIndices = rec8indices; |
728 |
ret.elementFactor = 6; |
729 |
// fall through |
730 |
case Hex20Face: |
731 |
case Rec8_Contact: |
732 |
case Hex20Face_Contact: |
733 |
//VTK_QUADRATIC_QUAD |
734 |
ret.elementSize = 3; |
735 |
ret.elementType = ZONETYPE_TRIANGLE; |
736 |
ret.reducedElementSize = 4; |
737 |
ret.reducedElementType = ZONETYPE_QUAD; |
738 |
break; |
739 |
|
740 |
case Tet10: |
741 |
//VTK_QUADRATIC_TETRA |
742 |
ret.multiCellIndices = tet10indices; |
743 |
ret.elementFactor = 8; |
744 |
ret.elementSize = ret.reducedElementSize = 4; |
745 |
ret.elementType = ret.reducedElementType = ZONETYPE_TET; |
746 |
break; |
747 |
|
748 |
case Hex20: |
749 |
//VTK_QUADRATIC_HEXAHEDRON |
750 |
ret.multiCellIndices = hex20indices; |
751 |
ret.elementFactor = 36; |
752 |
ret.elementSize = 3; |
753 |
ret.elementType = ZONETYPE_TRIANGLE; |
754 |
ret.reducedElementSize = 8; |
755 |
ret.reducedElementType = ZONETYPE_HEX; |
756 |
break; |
757 |
|
758 |
case Hex27: |
759 |
ret.multiCellIndices = hex27indices; |
760 |
ret.elementFactor = 8; |
761 |
// fall through |
762 |
case Hex8: |
763 |
ret.elementSize = ret.reducedElementSize = 8; |
764 |
ret.elementType = ret.reducedElementType = ZONETYPE_HEX; |
765 |
break; |
766 |
|
767 |
default: |
768 |
cerr << "WARNING: Unknown Finley Type " << typeId << endl; |
769 |
break; |
770 |
} |
771 |
return ret; |
772 |
} |
773 |
|
774 |
} // namespace EscriptReader |
775 |
|