/[escript]/trunk/dataexporter/src/ElementData.cpp
ViewVC logotype

Contents of /trunk/dataexporter/src/ElementData.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2849 - (show annotations)
Fri Jan 15 03:45:52 2010 UTC (9 years, 7 months ago) by caltinay
File size: 20587 byte(s)
Oops...

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

  ViewVC Help
Powered by ViewVC 1.1.26