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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2854 - (show annotations)
Mon Jan 18 01:54:04 2010 UTC (10 years, 9 months ago) by caltinay
File size: 20550 byte(s)
Fixed a compile error and removed unneeded member.

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

  ViewVC Help
Powered by ViewVC 1.1.26