/[escript]/trunk/tools/libescriptreader/src/escriptreader/ElementData.cpp
ViewVC logotype

Annotation of /trunk/tools/libescriptreader/src/escriptreader/ElementData.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2187 - (hide annotations)
Tue Dec 23 04:13:15 2008 UTC (12 years, 7 months ago) by caltinay
File size: 22907 byte(s)
Moved escriptreader related classes into own namespace.

1 caltinay 2183
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 caltinay 2187 namespace EscriptReader {
88    
89 caltinay 2183 //
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     return *(IntVec*)(NULL);
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     ret = DBPutUcdmesh(dbfile, siloMeshName.c_str(),
564     originalMesh->getNumDims(), NULL, &fullMesh->coords[0],
565     fullMesh->getNumNodes(), count, varName.c_str(),
566     /*"facelist"*/NULL, DB_FLOAT, NULL);
567    
568     // Point mesh is useful for debugging
569     //DBPutPointmesh(dbfile, "/pointmesh",
570     // originalMesh->getNumDims(), &fullMesh->coords[0],
571     // fullMesh->getNumNodes(), DB_FLOAT, NULL);
572    
573     if (ret != 0)
574     return false;
575    
576     // decide whether to additionally write out the reduced mesh
577     if (reducedCount > 0) {
578     if (siloPath == "/")
579     siloMeshName = siloPath + reducedMesh->getName();
580     else
581     siloMeshName = siloPath + string("/") + reducedMesh->getName();
582     arraylen = reducedCount * reducedNodesPerElement;
583     eltype = toSiloElementType(reducedType);
584     varName = string("Reduced") + name + string("_zones");
585     ret = DBPutZonelist2(dbfile, varName.c_str(), reducedCount,
586     originalMesh->getNumDims(), &reducedNodes[0], arraylen, 0, 0,
587     numReducedGhostElements, &eltype, &reducedNodesPerElement,
588     &reducedCount, 1, NULL);
589     if (ret == 0)
590     ret = DBPutUcdmesh(dbfile, siloMeshName.c_str(),
591     originalMesh->getNumDims(), NULL, &reducedMesh->coords[0],
592     reducedMesh->getNumNodes(), reducedCount, varName.c_str(),
593     NULL, DB_FLOAT, NULL);
594     if (ret != 0)
595     return false;
596     numCells = reducedCount;
597     } else {
598     numCells = count;
599     }
600     meshName = siloMeshName.c_str();
601    
602     // finally, write out the element-centered variables on the correct mesh
603     varName = name + string("_Color");
604     ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName,
605     (float*)&color[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL);
606     if (ret == 0) {
607     varName = name + string("_Id");
608     ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName,
609     (float*)&ID[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL);
610     }
611     if (ret == 0) {
612     varName = name + string("_Owner");
613     ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName,
614     (float*)&owner[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL);
615     }
616     if (ret == 0) {
617     varName = name + string("_Tag");
618     ret = DBPutUcdvar1(dbfile, varName.c_str(), meshName,
619     (float*)&tag[0], numCells, NULL, 0, DB_INT, DB_ZONECENT, NULL);
620     }
621    
622     DBSetDir(dbfile, "/");
623     return (ret == 0);
624    
625     #else // !HAVE_SILO
626     return false;
627     #endif
628     }
629    
630     //
631     //
632     //
633     FinleyElementInfo ElementData::getFinleyTypeInfo(ElementTypeId typeId)
634     {
635     FinleyElementInfo ret;
636     ret.multiCellIndices = NULL;
637     ret.elementFactor = 1;
638    
639     switch (typeId) {
640     case Point1_Contact://untested
641     case Line2Face_Contact://untested
642     case Line3Face_Contact://untested
643     case Line2Face://untested
644     case Line3Face://untested
645     case Point1://untested
646     cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
647     ret.elementSize = 1;
648     ret.elementType = ZONETYPE_POLYGON;
649     break;
650    
651     case Tri3Face_Contact://untested
652     case Tri3Face://untested
653     cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
654     case Line2_Contact:
655     case Rec4Face_Contact:
656     case Rec4Face:
657     case Line2:
658     ret.elementSize = ret.reducedElementSize = 2;
659     ret.elementType = ret.reducedElementType = ZONETYPE_BEAM;
660     break;
661    
662     case Line3:
663     ret.multiCellIndices = line3indices;
664     ret.elementFactor = 2;
665     // fall through
666     case Line3_Contact:
667     case Tri6Face_Contact://untested
668     case Rec8Face_Contact:
669     case Tri6Face://untested
670     case Rec8Face:
671     //VTK_QUADRATIC_EDGE
672     ret.elementSize = ret.reducedElementSize = 2;
673     ret.elementType = ret.reducedElementType = ZONETYPE_BEAM;
674     break;
675    
676     case Tet4Face_Contact://untested
677     case Tet4Face://untested
678     cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
679     case Tri3_Contact:
680     case Tri3:
681     ret.elementSize = ret.reducedElementSize = 3;
682     ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
683     break;
684    
685     case Rec4_Contact:
686     case Hex8Face_Contact:
687     case Hex8Face:
688     case Rec4:
689     ret.elementSize = ret.reducedElementSize = 4;
690     ret.elementType = ret.reducedElementType = ZONETYPE_QUAD;
691     break;
692    
693     case Rec9:
694     ret.multiCellIndices = rec9indices;
695     ret.elementFactor = 4;
696     // fall through
697     case Rec9_Contact:
698     ret.elementSize = ret.reducedElementSize = 4;
699     ret.elementType = ret.reducedElementType = ZONETYPE_QUAD;
700     break;
701    
702     case Tet4:
703     ret.elementSize = ret.reducedElementSize = 4;
704     ret.elementType = ret.reducedElementType = ZONETYPE_TET;
705     break;
706    
707     case Tri6:
708     ret.multiCellIndices = tri6indices;
709     ret.elementFactor = 4;
710     // fall through
711     case Tri6_Contact:
712     case Tet10Face_Contact://untested
713     case Tet10Face://untested
714     //VTK_QUADRATIC_TRIANGLE
715     ret.elementSize = ret.reducedElementSize = 3;
716     ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
717     break;
718    
719     case Rec8:
720     ret.multiCellIndices = rec8indices;
721     ret.elementFactor = 6;
722     // fall through
723     case Hex20Face:
724     case Rec8_Contact:
725     case Hex20Face_Contact:
726     //VTK_QUADRATIC_QUAD
727     ret.elementSize = 3;
728     ret.elementType = ZONETYPE_TRIANGLE;
729     ret.reducedElementSize = 4;
730     ret.reducedElementType = ZONETYPE_QUAD;
731     break;
732    
733     case Tet10:
734     //VTK_QUADRATIC_TETRA
735     ret.multiCellIndices = tet10indices;
736     ret.elementFactor = 8;
737     ret.elementSize = ret.reducedElementSize = 4;
738     ret.elementType = ret.reducedElementType = ZONETYPE_TET;
739     break;
740    
741     case Hex20:
742     //VTK_QUADRATIC_HEXAHEDRON
743     ret.multiCellIndices = hex20indices;
744     ret.elementFactor = 36;
745     ret.elementSize = 3;
746     ret.elementType = ZONETYPE_TRIANGLE;
747     ret.reducedElementSize = 8;
748     ret.reducedElementType = ZONETYPE_HEX;
749     break;
750    
751     case Hex27:
752     ret.multiCellIndices = hex27indices;
753     ret.elementFactor = 8;
754     // fall through
755     case Hex8:
756     ret.elementSize = ret.reducedElementSize = 8;
757     ret.elementType = ret.reducedElementType = ZONETYPE_HEX;
758     break;
759    
760     default:
761     cerr << "WARNING: Unknown Finley Type " << typeId << endl;
762     break;
763     }
764     return ret;
765     }
766    
767 caltinay 2187 } // namespace EscriptReader
768    

  ViewVC Help
Powered by ViewVC 1.1.26