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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2910 - (show annotations)
Wed Feb 3 03:22:31 2010 UTC (9 years, 6 months ago) by caltinay
File size: 28080 byte(s)
Added a libescriptreader.a target which does not depend on escript or finley
at runtime. This will be used for the VisIt plugin so is not built by default.

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

  ViewVC Help
Powered by ViewVC 1.1.26