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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2898 - (hide annotations)
Mon Feb 1 01:23:46 2010 UTC (9 years, 7 months ago) by caltinay
File size: 27811 byte(s)
Reverted the move of system_dep.h since it didn't solve the problem nicely.
Created a wrapper class for the needed functionality instead.

1 caltinay 2183
2     /*******************************************************
3     *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 caltinay 2183 * 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 caltinay 2810 #include <escriptexport/ElementData.h>
15     #include <escriptexport/NodeData.h>
16    
17 caltinay 2898 #include <finley/CppAdapter/MeshAdapter.h>
18 caltinay 2810
19 caltinay 2819 #include <iostream>
20    
21 caltinay 2810 #if USE_NETCDF
22 caltinay 2888 #include <netcdfcpp.h>
23 caltinay 2810 #endif
24    
25     #if USE_SILO
26 caltinay 2183 #include <silo.h>
27     #endif
28    
29     using namespace std;
30    
31     // The following arrays contain indices to convert unsupported element
32     // types into supported ones (lines, triangles, quads, hexahedrons)
33    
34     static const size_t line3indices[2*2] = {
35     0, 2,
36     2, 1
37     };
38     static const size_t tri6indices[4*3] = {
39     0, 3, 5,
40     3, 1, 4,
41     5, 4, 2,
42     3, 4, 5
43     };
44     static const size_t rec8indices[6*3] = {
45     0, 4, 7,
46     4, 1, 5,
47     5, 2, 6,
48     6, 3, 7,
49 caltinay 2810 7, 5, 6,
50     7, 4, 5
51 caltinay 2183 };
52     static const size_t rec9indices[4*4] = {
53     0, 4, 8, 7,
54     4, 1, 5, 8,
55     7, 8, 6, 3,
56     8, 5, 2, 6
57     };
58     static const size_t tet10indices[8*4] = {
59     6, 0, 7, 4,
60     2, 6, 9, 5,
61     9, 7, 3, 8,
62     4, 1, 8, 5,
63     6, 7, 9, 8,
64     5, 6, 9, 8,
65     6, 5, 4, 8,
66     7, 6, 4, 8
67     };
68     static const size_t hex20indices[36*3] = {
69     0, 8, 12, 8, 1, 13, 13, 5, 16,
70     16, 4, 12, 8, 13, 16, 8, 16, 12,
71     1, 9, 13, 9, 2, 14, 14, 6, 17,
72     17, 5, 13, 9, 14, 17, 9, 17, 13,
73     2, 10, 14, 10, 3, 15, 15, 7, 18,
74     18, 14, 6, 10, 15, 18, 10, 18, 14,
75     3, 11, 15, 11, 0, 12, 12, 4, 19,
76     19, 7, 15, 11, 12, 19, 11, 19, 15,
77     4, 16, 19, 16, 5, 17, 17, 6, 18,
78     18, 7, 19, 16, 17, 18, 16, 18, 19,
79     3, 10, 11, 10, 2, 9, 9, 1, 8,
80     8, 0, 11, 10, 9, 8, 10, 8, 11
81     };
82     static const size_t hex27indices[8*8] = {
83     0, 8, 20, 11, 12, 21, 26, 24,
84     8, 1, 9, 20, 21, 13, 22, 26,
85     11, 20, 10, 3, 24, 26, 23, 15,
86     20, 9, 2, 10, 26, 22, 14, 23,
87     12, 21, 26, 24, 4, 16, 25, 19,
88     21, 13, 22, 26, 16, 5, 17, 25,
89     24, 26, 23, 15, 19, 25, 18, 7,
90     26, 22, 14, 23, 25, 17, 6, 18
91     };
92    
93 caltinay 2810 namespace escriptexport {
94 caltinay 2187
95 caltinay 2183 //
96     // Constructor
97     //
98 caltinay 2810 ElementData::ElementData(const string& elementName, NodeData_ptr nodeData)
99 caltinay 2834 : originalMesh(nodeData), name(elementName), numElements(0),
100 caltinay 2854 numGhostElements(0), nodesPerElement(0),
101 caltinay 2880 type(ZONETYPE_UNKNOWN), finleyTypeId(NoRef), elementFactor(1)
102 caltinay 2183 {
103 caltinay 2880 nodeMesh.reset(new NodeData(name));
104 caltinay 2183 }
105    
106     //
107     // Copy constructor
108     //
109     ElementData::ElementData(const ElementData& e)
110     {
111     name = e.name;
112 caltinay 2810 numElements = e.numElements;
113 caltinay 2183 numGhostElements = e.numGhostElements;
114     type = e.type;
115 caltinay 2880 finleyTypeId = e.finleyTypeId;
116 caltinay 2183 nodesPerElement = e.nodesPerElement;
117 caltinay 2880 elementFactor = e.elementFactor;
118 caltinay 2834 originalMesh = e.originalMesh;
119 caltinay 2848 if (e.nodeMesh)
120     nodeMesh.reset(new NodeData(*e.nodeMesh));
121 caltinay 2880 else
122     nodeMesh.reset(new NodeData(name));
123 caltinay 2183
124     nodes = e.nodes;
125     ID = e.ID;
126     color = e.color;
127     tag = e.tag;
128     owner = e.owner;
129 caltinay 2834
130     if (e.reducedElements)
131     reducedElements = ElementData_ptr(new ElementData(*e.reducedElements));
132 caltinay 2183 }
133    
134     //
135     //
136 caltinay 2880 //
137 caltinay 2810 bool ElementData::initFromFinley(const Finley_ElementFile* finleyFile)
138     {
139     numElements = finleyFile->numElements;
140    
141     if (numElements > 0) {
142     nodesPerElement = finleyFile->numNodes;
143    
144     int* iPtr;
145    
146     iPtr = finleyFile->Nodes;
147     nodes.clear();
148     nodes.insert(nodes.end(), numElements*nodesPerElement, 0);
149     copy(iPtr, iPtr+numElements*nodesPerElement, nodes.begin());
150    
151     iPtr = finleyFile->Color;
152     color.clear();
153     color.insert(color.end(), numElements, 0);
154     copy(iPtr, iPtr+numElements, color.begin());
155    
156     iPtr = finleyFile->Id;
157     ID.clear();
158     ID.insert(ID.end(), numElements, 0);
159     copy(iPtr, iPtr+numElements, ID.begin());
160    
161     iPtr = finleyFile->Owner;
162     owner.clear();
163     owner.insert(owner.end(), numElements, 0);
164     copy(iPtr, iPtr+numElements, owner.begin());
165    
166     iPtr = finleyFile->Tag;
167     tag.clear();
168     tag.insert(tag.end(), numElements, 0);
169     copy(iPtr, iPtr+numElements, tag.begin());
170    
171 caltinay 2880 finleyTypeId = finleyFile->referenceElementSet->referenceElement
172     ->Type->TypeId;
173     FinleyElementInfo f = getFinleyTypeInfo(finleyTypeId);
174 caltinay 2810 type = f.elementType;
175 caltinay 2880 elementFactor = f.elementFactor;
176     if (elementFactor > 1 || f.reducedElementSize != nodesPerElement)
177 caltinay 2810 buildReducedElements(f);
178    
179 caltinay 2880 if (f.useQuadNodes) {
180     CoordArray quadNodes;
181     int numQuadNodes;
182     Finley_ShapeFunction* sf = finleyFile->referenceElementSet
183     ->referenceElement->Parametrization;
184     numQuadNodes = sf->numQuadNodes;
185     for (int i=0; i<f.quadDim; i++) {
186     double* srcPtr = sf->QuadNodes+i;
187     float* c = new float[numQuadNodes];
188     quadNodes.push_back(c);
189     for (int j=0; j<numQuadNodes; j++, srcPtr+=f.quadDim) {
190     *c++ = (float) *srcPtr;
191     }
192     }
193     quadMask = buildQuadMask(quadNodes, numQuadNodes);
194     for (int i=0; i<f.quadDim; i++)
195     delete[] quadNodes[i];
196     quadNodes.clear();
197    
198    
199     // now the reduced quadrature
200     sf = finleyFile->referenceElementSet
201     ->referenceElementReducedQuadrature->Parametrization;
202     numQuadNodes = sf->numQuadNodes;
203     for (int i=0; i<f.quadDim; i++) {
204     double* srcPtr = sf->QuadNodes+i;
205     float* c = new float[numQuadNodes];
206     quadNodes.push_back(c);
207     for (int j=0; j<numQuadNodes; j++, srcPtr+=f.quadDim) {
208     *c++ = (float) *srcPtr;
209     }
210     }
211     reducedQuadMask = buildQuadMask(quadNodes, numQuadNodes);
212     for (int i=0; i<f.quadDim; i++)
213     delete[] quadNodes[i];
214     quadNodes.clear();
215     }
216    
217 caltinay 2810 buildMeshes();
218     }
219     return true;
220     }
221    
222 caltinay 2183 //
223     // Reads element data from given NetCDF file
224     //
225     bool ElementData::readFromNc(NcFile* ncfile)
226     {
227 caltinay 2810 #if USE_NETCDF
228 caltinay 2183 string num_str("num_");
229     num_str += name;
230    
231 caltinay 2810 NcAtt* att = ncfile->get_att(num_str.c_str());
232     numElements = att->as_int(0);
233 caltinay 2183
234 caltinay 2848 // Only attempt to read further if there are any elements.
235 caltinay 2183 // Having no elements is not an error.
236 caltinay 2810 if (numElements > 0) {
237     att = ncfile->get_att((num_str + string("_numNodes")).c_str());
238     nodesPerElement = att->as_int(0);
239 caltinay 2183
240 caltinay 2810 nodes.insert(nodes.end(), numElements*nodesPerElement, 0);
241     NcVar* var = ncfile->get_var((name + string("_Nodes")).c_str());
242     var->get(&nodes[0], numElements, nodesPerElement);
243 caltinay 2183
244 caltinay 2810 color.insert(color.end(), numElements, 0);
245     var = ncfile->get_var((name + string("_Color")).c_str());
246     var->get(&color[0], numElements);
247 caltinay 2183
248 caltinay 2810 ID.insert(ID.end(), numElements, 0);
249     var = ncfile->get_var((name + string("_Id")).c_str());
250     var->get(&ID[0], numElements);
251 caltinay 2183
252 caltinay 2810 owner.insert(owner.end(), numElements, 0);
253     var = ncfile->get_var((name + string("_Owner")).c_str());
254     var->get(&owner[0], numElements);
255 caltinay 2183
256 caltinay 2810 tag.insert(tag.end(), numElements, 0);
257     var = ncfile->get_var((name + string("_Tag")).c_str());
258     var->get(&tag[0], numElements);
259 caltinay 2183
260 caltinay 2810 att = ncfile->get_att((name + string("_TypeId")).c_str());
261 caltinay 2880 finleyTypeId = (ElementTypeId)att->as_int(0);
262     FinleyElementInfo f = getFinleyTypeInfo(finleyTypeId);
263 caltinay 2810 type = f.elementType;
264 caltinay 2880 elementFactor = f.elementFactor;
265 caltinay 2810 if (f.elementFactor > 1 || f.reducedElementSize != nodesPerElement)
266     buildReducedElements(f);
267 caltinay 2183
268 caltinay 2880 if (f.useQuadNodes) {
269     att = ncfile->get_att("order");
270     int order = att->as_int(0);
271     att = ncfile->get_att("reduced_order");
272     int reduced_order = att->as_int(0);
273 caltinay 2898 finley::ReferenceElementSetWrapper wrapper(finleyTypeId, order,
274     reduced_order);
275     Finley_ReferenceElementSet* refElements = wrapper.getElementSet();
276 caltinay 2880
277     CoordArray quadNodes;
278     int numQuadNodes;
279     Finley_ShapeFunction* sf = refElements->referenceElement
280     ->Parametrization;
281     numQuadNodes = sf->numQuadNodes;
282     for (int i=0; i<f.quadDim; i++) {
283     double* srcPtr = sf->QuadNodes+i;
284     float* c = new float[numQuadNodes];
285     quadNodes.push_back(c);
286     for (int j=0; j<numQuadNodes; j++, srcPtr+=f.quadDim) {
287     *c++ = (float) *srcPtr;
288     }
289     }
290     quadMask = buildQuadMask(quadNodes, numQuadNodes);
291     for (int i=0; i<f.quadDim; i++)
292     delete[] quadNodes[i];
293     quadNodes.clear();
294    
295     // now the reduced quadrature
296     sf = refElements->referenceElementReducedQuadrature->Parametrization;
297     numQuadNodes = sf->numQuadNodes;
298     for (int i=0; i<f.quadDim; i++) {
299     double* srcPtr = sf->QuadNodes+i;
300     float* c = new float[numQuadNodes];
301     quadNodes.push_back(c);
302     for (int j=0; j<numQuadNodes; j++, srcPtr+=f.quadDim) {
303     *c++ = (float) *srcPtr;
304     }
305     }
306     reducedQuadMask = buildQuadMask(quadNodes, numQuadNodes);
307     for (int i=0; i<f.quadDim; i++)
308     delete[] quadNodes[i];
309     quadNodes.clear();
310     }
311    
312 caltinay 2810 buildMeshes();
313     }
314 caltinay 2183
315     return true;
316 caltinay 2810 #else // !USE_NETCDF
317     return false;
318     #endif
319 caltinay 2183 }
320    
321     //
322     //
323     //
324 caltinay 2880 StringVec ElementData::getMeshNames() const
325     {
326     StringVec res;
327     if (nodeMesh)
328     res.push_back(nodeMesh->getName());
329     if (reducedElements) {
330     StringVec rNames = reducedElements->getMeshNames();
331     if (rNames.size() > 0)
332     res.insert(res.end(), rNames.begin(), rNames.end());
333     }
334     return res;
335     }
336    
337     //
338     //
339     //
340     StringVec ElementData::getVarNames() const
341     {
342     StringVec res;
343     res.push_back(name + string("_Color"));
344     res.push_back(name + string("_Id"));
345     res.push_back(name + string("_Owner"));
346     res.push_back(name + string("_Tag"));
347     return res;
348     }
349    
350     //
351     //
352     //
353     const IntVec& ElementData::getVarDataByName(const string varName) const
354     {
355     if (varName == name+string("_Color"))
356     return color;
357     else if (varName == name+string("_Id"))
358     return ID;
359     else if (varName == name+string("_Owner")) {
360     return owner;
361     } else if (varName == name+string("_Tag"))
362     return tag;
363     else if (reducedElements)
364     return reducedElements->getVarDataByName(varName);
365     else
366     throw "Invalid variable name";
367     }
368    
369     //
370     //
371     //
372     const QuadMaskInfo& ElementData::getQuadMask(int functionSpace) const
373     {
374     if (functionSpace == FINLEY_REDUCED_ELEMENTS ||
375     functionSpace == FINLEY_REDUCED_FACE_ELEMENTS ||
376     functionSpace == FINLEY_REDUCED_CONTACT_ELEMENTS_1) {
377     return reducedQuadMask;
378     } else {
379     return quadMask;
380     }
381     }
382    
383     //
384     //
385     //
386 caltinay 2183 void ElementData::reorderArray(IntVec& v, const IntVec& idx,
387     int elementsPerIndex)
388     {
389     IntVec newArray(v.size());
390     IntVec::iterator arrIt = newArray.begin();
391     IntVec::const_iterator idxIt;
392     if (elementsPerIndex == 1) {
393     for (idxIt=idx.begin(); idxIt!=idx.end(); idxIt++) {
394     *arrIt++ = v[*idxIt];
395     }
396     } else {
397     for (idxIt=idx.begin(); idxIt!=idx.end(); idxIt++) {
398     int i = *idxIt;
399     copy(&v[i*elementsPerIndex], &v[(i+1)*elementsPerIndex], arrIt);
400     arrIt += elementsPerIndex;
401     }
402     }
403     v.swap(newArray);
404     }
405    
406     //
407     //
408     //
409     void ElementData::buildReducedElements(const FinleyElementInfo& f)
410     {
411 caltinay 2834 // create the node list for the new element type
412     IntVec reducedNodes(f.reducedElementSize*numElements, 0);
413    
414 caltinay 2183 IntVec::iterator reducedIt = reducedNodes.begin();
415     IntVec::const_iterator origIt;
416     for (origIt=nodes.begin(); origIt!=nodes.end();
417     origIt+=nodesPerElement)
418     {
419 caltinay 2834 copy(origIt, origIt+f.reducedElementSize, reducedIt);
420 caltinay 2183 reducedIt += f.reducedElementSize;
421     }
422    
423 caltinay 2834 if (f.elementFactor > 1) {
424     // replace each element by multiple smaller ones which will be the
425     // new 'full' elements, whereas the original ones are the reduced
426     // elements, e.g.:
427     //
428     // new reduced: new full:
429     // _________ _________
430     // | | | | |
431     // | | > |____|____|
432     // | | > | | |
433     // |_________| |____|____|
434     //
435    
436     // create the reduced elements which are basically a copy of the
437     // current elements
438     reducedElements = ElementData_ptr(new ElementData(
439     "Reduced"+name, originalMesh));
440     reducedElements->nodes = reducedNodes;
441     reducedElements->numElements = numElements;
442     reducedElements->type = f.reducedElementType;
443     reducedElements->nodesPerElement = f.reducedElementSize;
444     reducedElements->owner = owner;
445     reducedElements->color = color;
446     reducedElements->ID = ID;
447     reducedElements->tag = tag;
448    
449     // now update full element data
450 caltinay 2810 IntVec fullNodes(f.elementSize*f.elementFactor*numElements);
451 caltinay 2183 IntVec::iterator cellIt = fullNodes.begin();
452    
453     owner.clear();
454 caltinay 2834 color.clear();
455     ID.clear();
456     tag.clear();
457 caltinay 2810 for (int i=0; i < numElements; i++) {
458 caltinay 2834 owner.insert(owner.end(), f.elementFactor, reducedElements->owner[i]);
459     color.insert(color.end(), f.elementFactor, reducedElements->color[i]);
460     ID.insert(ID.end(), f.elementFactor, reducedElements->ID[i]);
461     tag.insert(tag.end(), f.elementFactor, reducedElements->tag[i]);
462 caltinay 2183 for (int j=0; j < f.elementFactor*f.elementSize; j++)
463 caltinay 2834 *cellIt++ = nodes[nodesPerElement*i+f.multiCellIndices[j]];
464 caltinay 2183 }
465 caltinay 2834
466 caltinay 2183 nodes.swap(fullNodes);
467     nodesPerElement = f.elementSize;
468 caltinay 2810 numElements *= f.elementFactor;
469 caltinay 2834
470 caltinay 2183 } else {
471 caltinay 2834 // we merely converted element types and don't need reduced elements
472     // so just replace node list and type
473 caltinay 2183 nodes.swap(reducedNodes);
474     nodesPerElement = f.reducedElementSize;
475     type = f.reducedElementType;
476     }
477     }
478    
479     //
480     //
481     //
482 caltinay 2834 IntVec ElementData::prepareGhostIndices(int ownIndex)
483 caltinay 2183 {
484 caltinay 2834 IntVec indexArray;
485 caltinay 2183 numGhostElements = 0;
486    
487     // move indices of "ghost zones" to the end to be able to reorder
488     // data accordingly
489 caltinay 2834 for (int i=0; i<numElements; i++) {
490 caltinay 2183 if (owner[i] == ownIndex)
491     indexArray.push_back(i);
492 caltinay 2834 }
493 caltinay 2183
494 caltinay 2834 for (int i=0; i<numElements; i++) {
495 caltinay 2183 if (owner[i] != ownIndex) {
496     numGhostElements++;
497     indexArray.push_back(i);
498     }
499 caltinay 2834 }
500     return indexArray;
501 caltinay 2183 }
502    
503     //
504     //
505     //
506 caltinay 2810 void ElementData::reorderGhostZones(int ownIndex)
507 caltinay 2183 {
508 caltinay 2834 IntVec indexArray = prepareGhostIndices(ownIndex);
509 caltinay 2183
510     // move "ghost data" to the end of the arrays
511     if (numGhostElements > 0) {
512     reorderArray(nodes, indexArray, nodesPerElement);
513     reorderArray(owner, indexArray, 1);
514 caltinay 2834 reorderArray(color, indexArray, 1);
515     reorderArray(ID, indexArray, 1);
516     reorderArray(tag, indexArray, 1);
517 caltinay 2183 }
518    
519 caltinay 2834 if (reducedElements)
520     reducedElements->reorderGhostZones(ownIndex);
521 caltinay 2183 }
522    
523     //
524     //
525     //
526 caltinay 2810 void ElementData::removeGhostZones(int ownIndex)
527 caltinay 2183 {
528 caltinay 2810 reorderGhostZones(ownIndex);
529    
530 caltinay 2183 if (numGhostElements > 0) {
531 caltinay 2810 numElements -= numGhostElements;
532     nodes.resize(numElements*nodesPerElement);
533     owner.resize(numElements);
534 caltinay 2834 color.resize(numElements);
535     ID.resize(numElements);
536     tag.resize(numElements);
537 caltinay 2183 numGhostElements = 0;
538     }
539    
540 caltinay 2834 if (reducedElements)
541     reducedElements->removeGhostZones(ownIndex);
542 caltinay 2183 }
543    
544     //
545     //
546     //
547     void ElementData::buildMeshes()
548     {
549 caltinay 2810 // build a new mesh containing only the required nodes
550     if (numElements > 0) {
551 caltinay 2880 if (nodeMesh && nodeMesh->getNumNodes() > 0) {
552 caltinay 2834 NodeData_ptr newMesh(new NodeData(nodeMesh, nodes, name));
553     nodeMesh.swap(newMesh);
554     } else {
555     nodeMesh.reset(new NodeData(originalMesh, nodes, name));
556     }
557 caltinay 2183 #ifdef _DEBUG
558 caltinay 2834 cout << nodeMesh->getName() << " has " << nodeMesh->getNumNodes()
559     << " nodes and " << numElements << " elements" << endl;
560 caltinay 2183 #endif
561 caltinay 2810 }
562 caltinay 2183
563 caltinay 2834 if (reducedElements)
564     reducedElements->buildMeshes();
565 caltinay 2183 }
566    
567 caltinay 2886 //
568     //
569     //
570     void ElementData::writeConnectivityVTK(ostream& os)
571     {
572     if (numElements > 0) {
573     const IntVec& gNI = nodeMesh->getGlobalNodeIndices();
574     IntVec::const_iterator it;
575     int count = 1;
576     for (it=nodes.begin(); it!=nodes.end(); it++, count++) {
577     os << gNI[*it];
578     if (count % nodesPerElement == 0)
579     os << endl;
580     else
581     os << " ";
582     }
583     }
584     }
585    
586 caltinay 2810 #if USE_SILO
587 caltinay 2183 //
588     //
589     //
590     inline int toSiloElementType(int type)
591     {
592     switch (type) {
593     case ZONETYPE_BEAM: return DB_ZONETYPE_BEAM;
594     case ZONETYPE_HEX: return DB_ZONETYPE_HEX;
595     case ZONETYPE_POLYGON: return DB_ZONETYPE_POLYGON;
596     case ZONETYPE_QUAD: return DB_ZONETYPE_QUAD;
597     case ZONETYPE_TET: return DB_ZONETYPE_TET;
598     case ZONETYPE_TRIANGLE: return DB_ZONETYPE_TRIANGLE;
599     }
600     return 0;
601     }
602     #endif
603    
604     //
605     //
606     //
607     bool ElementData::writeToSilo(DBfile* dbfile, const string& siloPath)
608     {
609 caltinay 2810 #if USE_SILO
610     if (numElements == 0)
611 caltinay 2183 return true;
612    
613     int ret;
614    
615     if (siloPath != "") {
616     ret = DBSetDir(dbfile, siloPath.c_str());
617     if (ret != 0)
618     return false;
619     }
620    
621     // write out the full mesh in any case
622 caltinay 2834 nodeMesh->setSiloPath(siloPath);
623     string siloMeshNameStr = nodeMesh->getFullSiloName();
624 caltinay 2810 const char* siloMeshName = siloMeshNameStr.c_str();
625     int arraylen = numElements * nodesPerElement;
626     int eltype = toSiloElementType(type);
627 caltinay 2183
628 caltinay 2834 string varName = name + string("_zones");
629 caltinay 2810 ret = DBPutZonelist2(dbfile, varName.c_str(), numElements,
630 caltinay 2834 nodeMesh->getNumDims(), &nodes[0], arraylen, 0, 0,
631 caltinay 2810 numGhostElements, &eltype, &nodesPerElement, &numElements, 1, NULL);
632 caltinay 2196 if (ret == 0) {
633 caltinay 2834 CoordArray& coordbase = const_cast<CoordArray&>(nodeMesh->getCoords());
634 caltinay 2810 ret = DBPutUcdmesh(dbfile, siloMeshName,
635 caltinay 2834 nodeMesh->getNumDims(), NULL, &coordbase[0],
636     nodeMesh->getNumNodes(), numElements, varName.c_str(),
637 caltinay 2183 /*"facelist"*/NULL, DB_FLOAT, NULL);
638 caltinay 2196 }
639 caltinay 2183
640     // Point mesh is useful for debugging
641 caltinay 2196 if (0) {
642 caltinay 2834 CoordArray& coordbase = const_cast<CoordArray&>(nodeMesh->getCoords());
643     DBPutPointmesh(dbfile, "/pointmesh", nodeMesh->getNumDims(),
644     &coordbase[0], nodeMesh->getNumNodes(), DB_FLOAT, NULL);
645 caltinay 2196 }
646 caltinay 2183
647     if (ret != 0)
648     return false;
649    
650 caltinay 2834 // write out the element-centered variables
651 caltinay 2183 varName = name + string("_Color");
652 caltinay 2810 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
653 caltinay 2834 (float*)&color[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
654     NULL);
655 caltinay 2183 if (ret == 0) {
656     varName = name + string("_Id");
657 caltinay 2810 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
658 caltinay 2834 (float*)&ID[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
659     NULL);
660 caltinay 2183 }
661     if (ret == 0) {
662     varName = name + string("_Owner");
663 caltinay 2810 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
664 caltinay 2834 (float*)&owner[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
665     NULL);
666 caltinay 2183 }
667     if (ret == 0) {
668     varName = name + string("_Tag");
669 caltinay 2810 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
670 caltinay 2834 (float*)&tag[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
671     NULL);
672 caltinay 2183 }
673 caltinay 2810
674 caltinay 2834 if (reducedElements) {
675     reducedElements->writeToSilo(dbfile, siloPath);
676     }
677    
678     // "Elements" is a special case
679 caltinay 2810 if (name == "Elements") {
680 caltinay 2834 nodeMesh->writeToSilo(dbfile);
681 caltinay 2810 }
682    
683 caltinay 2183 return (ret == 0);
684    
685 caltinay 2810 #else // !USE_SILO
686 caltinay 2183 return false;
687     #endif
688     }
689    
690     //
691     //
692     //
693     FinleyElementInfo ElementData::getFinleyTypeInfo(ElementTypeId typeId)
694     {
695     FinleyElementInfo ret;
696     ret.multiCellIndices = NULL;
697     ret.elementFactor = 1;
698 caltinay 2880 ret.useQuadNodes = false;
699     ret.quadDim = 0;
700 caltinay 2183
701     switch (typeId) {
702     case Point1_Contact://untested
703     case Line2Face_Contact://untested
704     case Line3Face_Contact://untested
705     case Line2Face://untested
706     case Line3Face://untested
707     case Point1://untested
708     cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
709     ret.elementSize = 1;
710     ret.elementType = ZONETYPE_POLYGON;
711     break;
712    
713     case Tri3Face_Contact://untested
714     case Tri3Face://untested
715     cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
716     case Line2_Contact:
717     case Rec4Face_Contact:
718     case Rec4Face:
719     case Line2:
720     ret.elementSize = ret.reducedElementSize = 2;
721     ret.elementType = ret.reducedElementType = ZONETYPE_BEAM;
722     break;
723    
724     case Line3:
725     ret.multiCellIndices = line3indices;
726     ret.elementFactor = 2;
727     // fall through
728     case Line3_Contact:
729     case Tri6Face_Contact://untested
730     case Rec8Face_Contact:
731     case Tri6Face://untested
732     case Rec8Face:
733     //VTK_QUADRATIC_EDGE
734     ret.elementSize = ret.reducedElementSize = 2;
735     ret.elementType = ret.reducedElementType = ZONETYPE_BEAM;
736     break;
737    
738     case Tet4Face_Contact://untested
739     case Tet4Face://untested
740     cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
741     case Tri3_Contact:
742     case Tri3:
743     ret.elementSize = ret.reducedElementSize = 3;
744     ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
745     break;
746    
747     case Rec4_Contact:
748     case Hex8Face_Contact:
749     case Hex8Face:
750     case Rec4:
751     ret.elementSize = ret.reducedElementSize = 4;
752     ret.elementType = ret.reducedElementType = ZONETYPE_QUAD;
753     break;
754    
755     case Rec9:
756 caltinay 2880 ret.useQuadNodes = true;
757     ret.quadDim = 2;
758 caltinay 2183 ret.multiCellIndices = rec9indices;
759     ret.elementFactor = 4;
760     // fall through
761     case Rec9_Contact:
762     ret.elementSize = ret.reducedElementSize = 4;
763     ret.elementType = ret.reducedElementType = ZONETYPE_QUAD;
764     break;
765    
766     case Tet4:
767     ret.elementSize = ret.reducedElementSize = 4;
768     ret.elementType = ret.reducedElementType = ZONETYPE_TET;
769     break;
770    
771     case Tri6:
772     ret.multiCellIndices = tri6indices;
773     ret.elementFactor = 4;
774     // fall through
775     case Tri6_Contact:
776     case Tet10Face_Contact://untested
777     case Tet10Face://untested
778     //VTK_QUADRATIC_TRIANGLE
779     ret.elementSize = ret.reducedElementSize = 3;
780     ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
781     break;
782    
783     case Rec8:
784     ret.multiCellIndices = rec8indices;
785     ret.elementFactor = 6;
786     // fall through
787     case Hex20Face:
788     case Rec8_Contact:
789     case Hex20Face_Contact:
790     //VTK_QUADRATIC_QUAD
791     ret.elementSize = 3;
792     ret.elementType = ZONETYPE_TRIANGLE;
793     ret.reducedElementSize = 4;
794     ret.reducedElementType = ZONETYPE_QUAD;
795     break;
796    
797     case Tet10:
798     //VTK_QUADRATIC_TETRA
799     ret.multiCellIndices = tet10indices;
800     ret.elementFactor = 8;
801     ret.elementSize = ret.reducedElementSize = 4;
802     ret.elementType = ret.reducedElementType = ZONETYPE_TET;
803     break;
804    
805     case Hex20:
806     //VTK_QUADRATIC_HEXAHEDRON
807     ret.multiCellIndices = hex20indices;
808     ret.elementFactor = 36;
809     ret.elementSize = 3;
810     ret.elementType = ZONETYPE_TRIANGLE;
811     ret.reducedElementSize = 8;
812     ret.reducedElementType = ZONETYPE_HEX;
813     break;
814    
815     case Hex27:
816 caltinay 2880 ret.useQuadNodes = true;
817     ret.quadDim = 3;
818 caltinay 2183 ret.multiCellIndices = hex27indices;
819     ret.elementFactor = 8;
820     // fall through
821     case Hex8:
822     ret.elementSize = ret.reducedElementSize = 8;
823     ret.elementType = ret.reducedElementType = ZONETYPE_HEX;
824     break;
825    
826     default:
827     cerr << "WARNING: Unknown Finley Type " << typeId << endl;
828     break;
829     }
830     return ret;
831     }
832    
833 caltinay 2880 //
834     //
835     //
836     inline bool inside1D(float x, float c, float r)
837     {
838     return (ABS(x-c) <= r);
839     }
840    
841     //
842     //
843     //
844     inline bool inside2D(float x, float y, float cx, float cy, float r)
845     {
846     return (inside1D(x, cx, r) && inside1D(y, cy, r));
847     }
848    
849     //
850     //
851     //
852     inline bool inside3D(float x, float y, float z,
853     float cx, float cy, float cz, float r)
854     {
855     return (inside2D(x, y, cx, cy, r) && inside1D(z, cz, r));
856     }
857    
858     //
859     //
860     //
861     QuadMaskInfo ElementData::buildQuadMask(const CoordArray& qnodes, int numQNodes)
862     {
863     QuadMaskInfo qmi;
864    
865     if (numQNodes == 0)
866     return qmi;
867    
868     if (finleyTypeId == Line3Macro) {
869     for (int i=0; i<elementFactor; i++) {
870     const float bounds[] = { 0.25, 0.75 };
871     IntVec m(numQNodes, 0);
872     int hits = 0;
873     for (size_t j=0; j<numQNodes; j++) {
874     if (inside1D(qnodes[0][j], bounds[i], .25)) {
875     m[j] = 1;
876     hits++;
877     }
878     }
879     qmi.mask.push_back(m);
880     if (hits == 0)
881     qmi.factor.push_back(1);
882     else
883     qmi.factor.push_back(hits);
884     }
885     } else if ((finleyTypeId == Rec9) || (finleyTypeId == Rec9Macro)) {
886     for (int i=0; i<elementFactor; i++) {
887     const float bounds[][2] = { { 0.25, 0.25 }, { 0.75, 0.25 },
888     { 0.25, 0.75 }, { 0.75, 0.75 } };
889     IntVec m(numQNodes, 0);
890     int hits = 0;
891     for (size_t j=0; j<numQNodes; j++) {
892     if (inside2D(qnodes[0][j], qnodes[1][j],
893     bounds[i][0], bounds[i][1], 0.25)) {
894     m[j] = 1;
895     hits++;
896     }
897     }
898     qmi.mask.push_back(m);
899     if (hits == 0)
900     qmi.factor.push_back(1);
901     else
902     qmi.factor.push_back(hits);
903     }
904     } else if ((finleyTypeId == Hex27) || (finleyTypeId == Hex27Macro) ){
905     for (int i=0; i<elementFactor; i++) {
906     const float bounds[][3] = {
907     { 0.25, 0.25, 0.25 }, { 0.75, 0.25, 0.25 },
908     { 0.25, 0.75, 0.25 }, { 0.75, 0.75, 0.25 },
909     { 0.25, 0.25, 0.75 }, { 0.75, 0.25, 0.75 },
910     { 0.25, 0.75, 0.75 }, { 0.75, 0.75, 0.75 } };
911     IntVec m(numQNodes, 0);
912     int hits = 0;
913     for (size_t j=0; j<numQNodes; j++) {
914     if (inside3D(qnodes[0][j], qnodes[1][j], qnodes[2][j],
915     bounds[i][0], bounds[i][1], bounds[i][2], 0.25)) {
916     m[j] = 1;
917     hits++;
918     }
919     }
920     qmi.mask.push_back(m);
921     if (hits == 0)
922     qmi.factor.push_back(1);
923     else
924     qmi.factor.push_back(hits);
925     }
926     }
927     return qmi;
928     }
929    
930 caltinay 2810 } // namespace escriptexport
931 caltinay 2187

  ViewVC Help
Powered by ViewVC 1.1.26