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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 4 months ago) by jfenwick
File size: 27407 byte(s)
Don't panic.
Updating copyright stamps

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

  ViewVC Help
Powered by ViewVC 1.1.26