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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2184 - (hide annotations)
Mon Dec 22 04:27:26 2008 UTC (12 years, 3 months ago) by caltinay
Original Path: trunk/tools/libescriptreader/src/escriptreader/ElementData.cpp
File size: 22977 byte(s)
Scons'ified escriptreader library and tools. Added checks and options for the
Silo library which is required for escript2silo and optional for the library.
Updated shake75 options file.

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     enum {
26     ZONETYPE_BEAM=1,
27     ZONETYPE_HEX,
28     ZONETYPE_POLYGON,
29     ZONETYPE_QUAD,
30     ZONETYPE_TET,
31 caltinay 2184 ZONETYPE_TRIANGLE
32 caltinay 2183 };
33    
34     // The following arrays contain indices to convert unsupported element
35     // types into supported ones (lines, triangles, quads, hexahedrons)
36    
37     static const size_t line3indices[2*2] = {
38     0, 2,
39     2, 1
40     };
41     static const size_t tri6indices[4*3] = {
42     0, 3, 5,
43     3, 1, 4,
44     5, 4, 2,
45     3, 4, 5
46     };
47     static const size_t rec8indices[6*3] = {
48     0, 4, 7,
49     4, 1, 5,
50     5, 2, 6,
51     6, 3, 7,
52     4, 6, 7,
53     4, 5, 6
54     };
55     static const size_t rec9indices[4*4] = {
56     0, 4, 8, 7,
57     4, 1, 5, 8,
58     7, 8, 6, 3,
59     8, 5, 2, 6
60     };
61     static const size_t tet10indices[8*4] = {
62     6, 0, 7, 4,
63     2, 6, 9, 5,
64     9, 7, 3, 8,
65     4, 1, 8, 5,
66     6, 7, 9, 8,
67     5, 6, 9, 8,
68     6, 5, 4, 8,
69     7, 6, 4, 8
70     };
71     static const size_t hex20indices[36*3] = {
72     0, 8, 12, 8, 1, 13, 13, 5, 16,
73     16, 4, 12, 8, 13, 16, 8, 16, 12,
74     1, 9, 13, 9, 2, 14, 14, 6, 17,
75     17, 5, 13, 9, 14, 17, 9, 17, 13,
76     2, 10, 14, 10, 3, 15, 15, 7, 18,
77     18, 14, 6, 10, 15, 18, 10, 18, 14,
78     3, 11, 15, 11, 0, 12, 12, 4, 19,
79     19, 7, 15, 11, 12, 19, 11, 19, 15,
80     4, 16, 19, 16, 5, 17, 17, 6, 18,
81     18, 7, 19, 16, 17, 18, 16, 18, 19,
82     3, 10, 11, 10, 2, 9, 9, 1, 8,
83     8, 0, 11, 10, 9, 8, 10, 8, 11
84     };
85     static const size_t hex27indices[8*8] = {
86     0, 8, 20, 11, 12, 21, 26, 24,
87     8, 1, 9, 20, 21, 13, 22, 26,
88     11, 20, 10, 3, 24, 26, 23, 15,
89     20, 9, 2, 10, 26, 22, 14, 23,
90     12, 21, 26, 24, 4, 16, 25, 19,
91     21, 13, 22, 26, 16, 5, 17, 25,
92     24, 26, 23, 15, 19, 25, 18, 7,
93     26, 22, 14, 23, 25, 17, 6, 18
94     };
95    
96     //
97     // Constructor
98     //
99     ElementData::ElementData(const string& elementName, const Mesh* mainMesh)
100     : name(elementName), count(0), reducedCount(0), numGhostElements(0),
101     numReducedGhostElements(0), fullMesh(NULL), reducedMesh(NULL),
102     originalMesh(mainMesh), fullMeshIsOriginalMesh(false)
103     {
104     }
105    
106     //
107     // Copy constructor
108     //
109     ElementData::ElementData(const ElementData& e)
110     {
111     name = e.name;
112     count = e.count;
113     reducedCount = e.reducedCount;
114     numGhostElements = e.numGhostElements;
115     numReducedGhostElements = e.numReducedGhostElements;
116     numDims = e.numDims;
117     type = e.type;
118     reducedType = e.reducedType;
119     nodesPerElement = e.nodesPerElement;
120     reducedNodesPerElement = e.reducedNodesPerElement;
121    
122     originalMesh = e.originalMesh;
123     fullMeshIsOriginalMesh = e.fullMeshIsOriginalMesh;
124     if (fullMeshIsOriginalMesh)
125     fullMesh = const_cast<Mesh*>(originalMesh);
126     else if (e.fullMesh)
127     fullMesh = new Mesh(*e.fullMesh);
128     else
129     fullMesh = NULL;
130    
131     if (e.reducedMesh)
132     reducedMesh = new Mesh(*e.reducedMesh);
133     else
134     reducedMesh = NULL;
135    
136     nodes = e.nodes;
137     reducedNodes = e.reducedNodes;
138     ID = e.ID;
139     color = e.color;
140     tag = e.tag;
141     owner = e.owner;
142     reducedOwner = e.reducedOwner;
143     indexArray = e.indexArray;
144     reducedIndexArray = e.reducedIndexArray;
145     ID2idx = e.ID2idx;
146     }
147    
148     //
149     // Destructor
150     //
151     ElementData::~ElementData()
152     {
153     if (reducedMesh)
154     delete reducedMesh;
155     if (fullMesh && !fullMeshIsOriginalMesh)
156     delete fullMesh;
157     }
158    
159     StringVec ElementData::getMeshNames() const
160     {
161     StringVec res;
162     if (fullMesh && !fullMeshIsOriginalMesh)
163     res.push_back(fullMesh->getName());
164     if (reducedMesh)
165     res.push_back(reducedMesh->getName());
166     return res;
167     }
168    
169     StringVec ElementData::getVarNames() const
170     {
171     StringVec res;
172     if (count > 0) {
173     res.push_back(name + string("_Color"));
174     res.push_back(name + string("_Id"));
175     res.push_back(name + string("_Owner"));
176     res.push_back(name + string("_Tag"));
177     }
178     return res;
179     }
180    
181     const IntVec& ElementData::getVarDataByName(const string varName) const
182     {
183     if (varName == name+string("_Color"))
184     return color;
185     else if (varName == name+string("_Id"))
186     return ID;
187     else if (varName == name+string("_Owner")) {
188     if (reducedCount > 0)
189     return reducedOwner;
190     else
191     return owner;
192     } else if (varName == name+string("_Tag"))
193     return tag;
194     else
195     return *(IntVec*)(NULL);
196     }
197    
198     //
199     // Reads element data from given NetCDF file
200     //
201     bool ElementData::readFromNc(NcFile* ncfile)
202     {
203     NcAtt* att;
204     NcVar* var;
205    
206     string num_str("num_");
207     num_str += name;
208    
209     att = ncfile->get_att(num_str.c_str());
210     count = att->as_int(0);
211    
212     // Only attempt to read data if there are any elements.
213     // Having no elements is not an error.
214     if (count == 0)
215     return true;
216    
217     att = ncfile->get_att((num_str + string("_numNodes")).c_str());
218     nodesPerElement = att->as_int(0);
219    
220     nodes.insert(nodes.end(), count*nodesPerElement, 0);
221     var = ncfile->get_var((name + string("_Nodes")).c_str());
222     var->get(&nodes[0], count, nodesPerElement);
223    
224     color.insert(color.end(), count, 0);
225     var = ncfile->get_var((name + string("_Color")).c_str());
226     var->get(&color[0], count);
227    
228     ID.insert(ID.end(), count, 0);
229     var = ncfile->get_var((name + string("_Id")).c_str());
230     var->get(&ID[0], count);
231    
232     owner.insert(owner.end(), count, 0);
233     var = ncfile->get_var((name + string("_Owner")).c_str());
234     var->get(&owner[0], count);
235    
236     tag.insert(tag.end(), count, 0);
237     var = ncfile->get_var((name + string("_Tag")).c_str());
238     var->get(&tag[0], count);
239    
240     att = ncfile->get_att((name + string("_TypeId")).c_str());
241     FinleyElementInfo f = getFinleyTypeInfo((ElementTypeId)att->as_int(0));
242     type = f.elementType;
243     reducedType = f.reducedElementType;
244    
245     // build elementID->index map
246     buildIndexMap();
247    
248     if (f.elementFactor > 1 || f.reducedElementSize != nodesPerElement)
249     buildReducedElements(f);
250    
251     buildMeshes();
252    
253     return true;
254     }
255    
256     //
257     //
258     //
259     void ElementData::reorderArray(IntVec& v, const IntVec& idx,
260     int elementsPerIndex)
261     {
262     IntVec newArray(v.size());
263     IntVec::iterator arrIt = newArray.begin();
264     IntVec::const_iterator idxIt;
265     if (elementsPerIndex == 1) {
266     for (idxIt=idx.begin(); idxIt!=idx.end(); idxIt++) {
267     *arrIt++ = v[*idxIt];
268     }
269     } else {
270     for (idxIt=idx.begin(); idxIt!=idx.end(); idxIt++) {
271     int i = *idxIt;
272     copy(&v[i*elementsPerIndex], &v[(i+1)*elementsPerIndex], arrIt);
273     arrIt += elementsPerIndex;
274     }
275     }
276     v.swap(newArray);
277     }
278    
279     //
280     //
281     //
282     void ElementData::buildReducedElements(const FinleyElementInfo& f)
283     {
284     reducedNodes.clear();
285     reducedNodes.insert(reducedNodes.end(),
286     f.reducedElementSize*count, 0);
287     IntVec::iterator reducedIt = reducedNodes.begin();
288     IntVec::const_iterator origIt;
289     for (origIt=nodes.begin(); origIt!=nodes.end();
290     origIt+=nodesPerElement)
291     {
292     std::copy(origIt, origIt+f.reducedElementSize, reducedIt);
293     reducedIt += f.reducedElementSize;
294     }
295    
296     // Remove comment to save disk space - we don't really need the full
297     // elements except for the main mesh
298     if (f.elementFactor > 1 /*&& name == "Elements"*/) {
299     // replace each element by multiple smaller ones
300     IntVec fullNodes(f.elementSize*f.elementFactor*count);
301     IntVec::iterator cellIt = fullNodes.begin();
302    
303     // copy owner data
304     owner.swap(reducedOwner);
305     owner.clear();
306     for (int i=0; i < count; i++) {
307     owner.insert(owner.end(), f.elementFactor, reducedOwner[i]);
308     for (int j=0; j < f.elementFactor*f.elementSize; j++)
309     *cellIt++ = nodes[
310     nodesPerElement*i+f.multiCellIndices[j]];
311     }
312     nodes.swap(fullNodes);
313     reducedCount = count;
314     reducedNodesPerElement = f.reducedElementSize;
315     nodesPerElement = f.elementSize;
316     count *= f.elementFactor;
317     } else {
318     // we only keep the reduced elements but treat them as regular
319     // ones, so replace the data accordingly
320     nodes.swap(reducedNodes);
321     reducedNodes.clear();
322     nodesPerElement = f.reducedElementSize;
323     type = f.reducedElementType;
324     }
325     }
326    
327     //
328     //
329     //
330     void ElementData::prepareGhostIndices(int ownIndex)
331     {
332     indexArray.clear();
333     reducedIndexArray.clear();
334     numGhostElements = 0;
335     numReducedGhostElements = 0;
336    
337     // move indices of "ghost zones" to the end to be able to reorder
338     // data accordingly
339     for (int i=0; i<count; i++)
340     if (owner[i] == ownIndex)
341     indexArray.push_back(i);
342    
343     for (int i=0; i<count; i++)
344     if (owner[i] != ownIndex) {
345     numGhostElements++;
346     indexArray.push_back(i);
347     }
348    
349     for (int i=0; i<reducedCount; i++)
350     if (reducedOwner[i] == ownIndex)
351     reducedIndexArray.push_back(i);
352    
353     for (int i=0; i<reducedCount; i++)
354     if (reducedOwner[i] != ownIndex) {
355     numReducedGhostElements++;
356     reducedIndexArray.push_back(i);
357     }
358     }
359    
360     //
361     //
362     //
363     void ElementData::handleGhostZones(int ownIndex)
364     {
365     prepareGhostIndices(ownIndex);
366    
367     // move "ghost data" to the end of the arrays
368     if (numGhostElements > 0) {
369     reorderArray(nodes, indexArray, nodesPerElement);
370     reorderArray(owner, indexArray, 1);
371     if (reducedCount == 0) {
372     reorderArray(color, indexArray, 1);
373     reorderArray(ID, indexArray, 1);
374     reorderArray(tag, indexArray, 1);
375     }
376     }
377    
378     if (numReducedGhostElements > 0) {
379     reorderArray(reducedNodes, reducedIndexArray, reducedNodesPerElement);
380     reorderArray(reducedOwner, reducedIndexArray, 1);
381     reorderArray(color, reducedIndexArray, 1);
382     reorderArray(ID, reducedIndexArray, 1);
383     reorderArray(tag, reducedIndexArray, 1);
384     }
385     }
386    
387     //
388     //
389     //
390     void ElementData::removeGhostZones()
391     {
392     if (numGhostElements > 0) {
393     count -= numGhostElements;
394     nodes.resize(count*nodesPerElement);
395     owner.resize(count);
396     if (reducedCount == 0) {
397     color.resize(count);
398     ID.resize(count);
399     tag.resize(count);
400     }
401     numGhostElements = 0;
402     }
403    
404     if (numReducedGhostElements > 0) {
405     reducedCount -= numReducedGhostElements;
406     reducedNodes.resize(reducedCount*reducedNodesPerElement);
407     reducedOwner.resize(reducedCount);
408     color.resize(reducedCount);
409     ID.resize(reducedCount);
410     tag.resize(reducedCount);
411     numReducedGhostElements = 0;
412     }
413     }
414    
415     //
416     //
417     //
418     void ElementData::buildMeshes()
419     {
420     if (count == 0)
421     return;
422    
423     // use existing original mesh for Elements but build a new mesh
424     // containing only the required nodes for other types
425     if (name == "Elements") {
426     fullMesh = const_cast<Mesh*>(originalMesh);
427     fullMeshIsOriginalMesh = true;
428     } else {
429     // first: build a map of required IDs while translating
430     // the original node IDs at the same time
431     IntVec::iterator it;
432     IndexMap nodeID2idx;
433     IntVec nodeIDs;
434     size_t newIdx = 0;
435    
436     for (it = nodes.begin(); it != nodes.end(); it++) {
437     IndexMap::iterator res = nodeID2idx.find(*it);
438     if (res == nodeID2idx.end()) {
439     nodeIDs.push_back(*it);
440     nodeID2idx[*it] = newIdx;
441     *it = newIdx++;
442     } else {
443     *it = res->second;
444     }
445     }
446    
447     // second: use map to fill coordinates for the nodes
448     // newIdx contains the number of nodes
449     CoordArray coords;
450     const CoordArray& origCoords = originalMesh->getCoords();
451     for (int dim=0; dim < originalMesh->getNumDims(); dim++) {
452     float* c = new float[newIdx];
453     coords.push_back(c);
454     IndexMap::const_iterator mIt;
455     for (mIt = nodeID2idx.begin(); mIt != nodeID2idx.end(); mIt++)
456     c[mIt->second] = origCoords[dim][mIt->first];
457     }
458    
459     if (fullMesh)
460     delete fullMesh;
461    
462     fullMesh = new Mesh(coords, originalMesh->getNumDims(), newIdx);
463     fullMesh->setName(name);
464     fullMesh->setIndexMap(nodeID2idx);
465     fullMesh->setNodeIDs(nodeIDs);
466     }
467    
468     #ifdef _DEBUG
469     cout << fullMesh->getName() << " has " << fullMesh->getNumNodes()
470     << " nodes, " << count << " elements" << endl;
471     #endif
472    
473     // build a reduced mesh if necessary
474     if (reducedCount > 0) {
475     IntVec::iterator it;
476     IndexMap nodeID2idx;
477     IntVec reducedNodeIDs;
478     size_t newIdx = 0;
479    
480     for (it = reducedNodes.begin(); it != reducedNodes.end(); it++) {
481     int id = originalMesh->getNodeIDs()[*it];
482     IndexMap::iterator res = nodeID2idx.find(id);
483     if (res == nodeID2idx.end()) {
484     reducedNodeIDs.push_back(id);
485     nodeID2idx[id] = newIdx;
486     *it = newIdx++;
487     } else {
488     *it = res->second;
489     }
490     }
491    
492     CoordArray coords;
493     const CoordArray& origCoords = originalMesh->getCoords();
494     for (int dim=0; dim < originalMesh->getNumDims(); dim++) {
495     float* c = new float[newIdx];
496     coords.push_back(c);
497     IndexMap::const_iterator mIt;
498     for (mIt = nodeID2idx.begin(); mIt != nodeID2idx.end(); mIt++) {
499     IndexMap::const_iterator idx;
500     idx = originalMesh->getIndexMap().find(mIt->first);
501     c[mIt->second] = origCoords[dim][idx->second];
502     }
503     }
504     if (reducedMesh)
505     delete reducedMesh;
506    
507     reducedMesh = new Mesh(coords, originalMesh->getNumDims(), newIdx);
508     reducedMesh->setName(string("Reduced") + name);
509     reducedMesh->setIndexMap(nodeID2idx);
510     reducedMesh->setNodeIDs(reducedNodeIDs);
511    
512     #ifdef _DEBUG
513     cout << reducedMesh->getName() << " has " << newIdx
514     << " nodes, " << reducedCount << " elements" << endl;
515     #endif
516     }
517     }
518    
519     #if HAVE_SILO
520     //
521     //
522     //
523     inline int toSiloElementType(int type)
524     {
525     switch (type) {
526     case ZONETYPE_BEAM: return DB_ZONETYPE_BEAM;
527     case ZONETYPE_HEX: return DB_ZONETYPE_HEX;
528     case ZONETYPE_POLYGON: return DB_ZONETYPE_POLYGON;
529     case ZONETYPE_QUAD: return DB_ZONETYPE_QUAD;
530     case ZONETYPE_TET: return DB_ZONETYPE_TET;
531     case ZONETYPE_TRIANGLE: return DB_ZONETYPE_TRIANGLE;
532     }
533     return 0;
534     }
535     #endif
536    
537     //
538     //
539     //
540     bool ElementData::writeToSilo(DBfile* dbfile, const string& siloPath)
541     {
542     #if HAVE_SILO
543     if (count == 0)
544     return true;
545    
546     const char* meshName;
547     int numCells;
548     string varName, siloMeshName;
549     int ret;
550    
551     if (siloPath != "") {
552     ret = DBSetDir(dbfile, siloPath.c_str());
553     if (ret != 0)
554     return false;
555     }
556    
557     // write out the full mesh in any case
558     if (siloPath == "/")
559     siloMeshName = siloPath + fullMesh->getName();
560     else
561     siloMeshName = siloPath + string("/") + fullMesh->getName();
562    
563     int arraylen = count * nodesPerElement;
564     int eltype = toSiloElementType(type);
565     varName = name + string("_zones");
566     ret = DBPutZonelist2(dbfile, varName.c_str(), count,
567     originalMesh->getNumDims(), &nodes[0], arraylen, 0, 0,
568     numGhostElements, &eltype, &nodesPerElement, &count, 1, NULL);
569     if (ret == 0)
570     ret = DBPutUcdmesh(dbfile, siloMeshName.c_str(),
571     originalMesh->getNumDims(), NULL, &fullMesh->coords[0],
572     fullMesh->getNumNodes(), count, varName.c_str(),
573     /*"facelist"*/NULL, DB_FLOAT, NULL);
574    
575     // Point mesh is useful for debugging
576     //DBPutPointmesh(dbfile, "/pointmesh",
577     // originalMesh->getNumDims(), &fullMesh->coords[0],
578     // fullMesh->getNumNodes(), DB_FLOAT, NULL);
579    
580     if (ret != 0)
581     return false;
582    
583     // decide whether to additionally write out the reduced mesh
584     if (reducedCount > 0) {
585     if (siloPath == "/")
586     siloMeshName = siloPath + reducedMesh->getName();
587     else
588     siloMeshName = siloPath + string("/") + reducedMesh->getName();
589     arraylen = reducedCount * reducedNodesPerElement;
590     eltype = toSiloElementType(reducedType);
591     varName = string("Reduced") + name + string("_zones");
592     ret = DBPutZonelist2(dbfile, varName.c_str(), reducedCount,
593     originalMesh->getNumDims(), &reducedNodes[0], arraylen, 0, 0,
594     numReducedGhostElements, &eltype, &reducedNodesPerElement,
595     &reducedCount, 1, NULL);
596     if (ret == 0)
597     ret = DBPutUcdmesh(dbfile, siloMeshName.c_str(),
598     originalMesh->getNumDims(), NULL, &reducedMesh->coords[0],
599     reducedMesh->getNumNodes(), reducedCount, varName.c_str(),
600     NULL, DB_FLOAT, NULL);
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    

  ViewVC Help
Powered by ViewVC 1.1.26