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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2898 - (show annotations)
Mon Feb 1 01:23:46 2010 UTC (9 years 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
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 #include <finley/CppAdapter/MeshAdapter.h>
18
19 #include <iostream>
20
21 #if USE_NETCDF
22 #include <netcdfcpp.h>
23 #endif
24
25 #if USE_SILO
26 #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 7, 5, 6,
50 7, 4, 5
51 };
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 namespace escriptexport {
94
95 //
96 // Constructor
97 //
98 ElementData::ElementData(const string& elementName, NodeData_ptr nodeData)
99 : originalMesh(nodeData), name(elementName), numElements(0),
100 numGhostElements(0), nodesPerElement(0),
101 type(ZONETYPE_UNKNOWN), finleyTypeId(NoRef), elementFactor(1)
102 {
103 nodeMesh.reset(new NodeData(name));
104 }
105
106 //
107 // Copy constructor
108 //
109 ElementData::ElementData(const ElementData& e)
110 {
111 name = e.name;
112 numElements = e.numElements;
113 numGhostElements = e.numGhostElements;
114 type = e.type;
115 finleyTypeId = e.finleyTypeId;
116 nodesPerElement = e.nodesPerElement;
117 elementFactor = e.elementFactor;
118 originalMesh = e.originalMesh;
119 if (e.nodeMesh)
120 nodeMesh.reset(new NodeData(*e.nodeMesh));
121 else
122 nodeMesh.reset(new NodeData(name));
123
124 nodes = e.nodes;
125 ID = e.ID;
126 color = e.color;
127 tag = e.tag;
128 owner = e.owner;
129
130 if (e.reducedElements)
131 reducedElements = ElementData_ptr(new ElementData(*e.reducedElements));
132 }
133
134 //
135 //
136 //
137 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 finleyTypeId = finleyFile->referenceElementSet->referenceElement
172 ->Type->TypeId;
173 FinleyElementInfo f = getFinleyTypeInfo(finleyTypeId);
174 type = f.elementType;
175 elementFactor = f.elementFactor;
176 if (elementFactor > 1 || f.reducedElementSize != nodesPerElement)
177 buildReducedElements(f);
178
179 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 buildMeshes();
218 }
219 return true;
220 }
221
222 //
223 // Reads element data from given NetCDF file
224 //
225 bool ElementData::readFromNc(NcFile* ncfile)
226 {
227 #if USE_NETCDF
228 string num_str("num_");
229 num_str += name;
230
231 NcAtt* att = ncfile->get_att(num_str.c_str());
232 numElements = att->as_int(0);
233
234 // Only attempt to read further if there are any elements.
235 // Having no elements is not an error.
236 if (numElements > 0) {
237 att = ncfile->get_att((num_str + string("_numNodes")).c_str());
238 nodesPerElement = att->as_int(0);
239
240 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
244 color.insert(color.end(), numElements, 0);
245 var = ncfile->get_var((name + string("_Color")).c_str());
246 var->get(&color[0], numElements);
247
248 ID.insert(ID.end(), numElements, 0);
249 var = ncfile->get_var((name + string("_Id")).c_str());
250 var->get(&ID[0], numElements);
251
252 owner.insert(owner.end(), numElements, 0);
253 var = ncfile->get_var((name + string("_Owner")).c_str());
254 var->get(&owner[0], numElements);
255
256 tag.insert(tag.end(), numElements, 0);
257 var = ncfile->get_var((name + string("_Tag")).c_str());
258 var->get(&tag[0], numElements);
259
260 att = ncfile->get_att((name + string("_TypeId")).c_str());
261 finleyTypeId = (ElementTypeId)att->as_int(0);
262 FinleyElementInfo f = getFinleyTypeInfo(finleyTypeId);
263 type = f.elementType;
264 elementFactor = f.elementFactor;
265 if (f.elementFactor > 1 || f.reducedElementSize != nodesPerElement)
266 buildReducedElements(f);
267
268 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 finley::ReferenceElementSetWrapper wrapper(finleyTypeId, order,
274 reduced_order);
275 Finley_ReferenceElementSet* refElements = wrapper.getElementSet();
276
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 buildMeshes();
313 }
314
315 return true;
316 #else // !USE_NETCDF
317 return false;
318 #endif
319 }
320
321 //
322 //
323 //
324 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 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 // create the node list for the new element type
412 IntVec reducedNodes(f.reducedElementSize*numElements, 0);
413
414 IntVec::iterator reducedIt = reducedNodes.begin();
415 IntVec::const_iterator origIt;
416 for (origIt=nodes.begin(); origIt!=nodes.end();
417 origIt+=nodesPerElement)
418 {
419 copy(origIt, origIt+f.reducedElementSize, reducedIt);
420 reducedIt += f.reducedElementSize;
421 }
422
423 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 IntVec fullNodes(f.elementSize*f.elementFactor*numElements);
451 IntVec::iterator cellIt = fullNodes.begin();
452
453 owner.clear();
454 color.clear();
455 ID.clear();
456 tag.clear();
457 for (int i=0; i < numElements; i++) {
458 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 for (int j=0; j < f.elementFactor*f.elementSize; j++)
463 *cellIt++ = nodes[nodesPerElement*i+f.multiCellIndices[j]];
464 }
465
466 nodes.swap(fullNodes);
467 nodesPerElement = f.elementSize;
468 numElements *= f.elementFactor;
469
470 } else {
471 // we merely converted element types and don't need reduced elements
472 // so just replace node list and type
473 nodes.swap(reducedNodes);
474 nodesPerElement = f.reducedElementSize;
475 type = f.reducedElementType;
476 }
477 }
478
479 //
480 //
481 //
482 IntVec ElementData::prepareGhostIndices(int ownIndex)
483 {
484 IntVec indexArray;
485 numGhostElements = 0;
486
487 // move indices of "ghost zones" to the end to be able to reorder
488 // data accordingly
489 for (int i=0; i<numElements; i++) {
490 if (owner[i] == ownIndex)
491 indexArray.push_back(i);
492 }
493
494 for (int i=0; i<numElements; i++) {
495 if (owner[i] != ownIndex) {
496 numGhostElements++;
497 indexArray.push_back(i);
498 }
499 }
500 return indexArray;
501 }
502
503 //
504 //
505 //
506 void ElementData::reorderGhostZones(int ownIndex)
507 {
508 IntVec indexArray = prepareGhostIndices(ownIndex);
509
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 reorderArray(color, indexArray, 1);
515 reorderArray(ID, indexArray, 1);
516 reorderArray(tag, indexArray, 1);
517 }
518
519 if (reducedElements)
520 reducedElements->reorderGhostZones(ownIndex);
521 }
522
523 //
524 //
525 //
526 void ElementData::removeGhostZones(int ownIndex)
527 {
528 reorderGhostZones(ownIndex);
529
530 if (numGhostElements > 0) {
531 numElements -= numGhostElements;
532 nodes.resize(numElements*nodesPerElement);
533 owner.resize(numElements);
534 color.resize(numElements);
535 ID.resize(numElements);
536 tag.resize(numElements);
537 numGhostElements = 0;
538 }
539
540 if (reducedElements)
541 reducedElements->removeGhostZones(ownIndex);
542 }
543
544 //
545 //
546 //
547 void ElementData::buildMeshes()
548 {
549 // build a new mesh containing only the required nodes
550 if (numElements > 0) {
551 if (nodeMesh && nodeMesh->getNumNodes() > 0) {
552 NodeData_ptr newMesh(new NodeData(nodeMesh, nodes, name));
553 nodeMesh.swap(newMesh);
554 } else {
555 nodeMesh.reset(new NodeData(originalMesh, nodes, name));
556 }
557 #ifdef _DEBUG
558 cout << nodeMesh->getName() << " has " << nodeMesh->getNumNodes()
559 << " nodes and " << numElements << " elements" << endl;
560 #endif
561 }
562
563 if (reducedElements)
564 reducedElements->buildMeshes();
565 }
566
567 //
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 #if USE_SILO
587 //
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 #if USE_SILO
610 if (numElements == 0)
611 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 nodeMesh->setSiloPath(siloPath);
623 string siloMeshNameStr = nodeMesh->getFullSiloName();
624 const char* siloMeshName = siloMeshNameStr.c_str();
625 int arraylen = numElements * nodesPerElement;
626 int eltype = toSiloElementType(type);
627
628 string varName = name + string("_zones");
629 ret = DBPutZonelist2(dbfile, varName.c_str(), numElements,
630 nodeMesh->getNumDims(), &nodes[0], arraylen, 0, 0,
631 numGhostElements, &eltype, &nodesPerElement, &numElements, 1, NULL);
632 if (ret == 0) {
633 CoordArray& coordbase = const_cast<CoordArray&>(nodeMesh->getCoords());
634 ret = DBPutUcdmesh(dbfile, siloMeshName,
635 nodeMesh->getNumDims(), NULL, &coordbase[0],
636 nodeMesh->getNumNodes(), numElements, varName.c_str(),
637 /*"facelist"*/NULL, DB_FLOAT, NULL);
638 }
639
640 // Point mesh is useful for debugging
641 if (0) {
642 CoordArray& coordbase = const_cast<CoordArray&>(nodeMesh->getCoords());
643 DBPutPointmesh(dbfile, "/pointmesh", nodeMesh->getNumDims(),
644 &coordbase[0], nodeMesh->getNumNodes(), DB_FLOAT, NULL);
645 }
646
647 if (ret != 0)
648 return false;
649
650 // write out the element-centered variables
651 varName = name + string("_Color");
652 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
653 (float*)&color[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
654 NULL);
655 if (ret == 0) {
656 varName = name + string("_Id");
657 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
658 (float*)&ID[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
659 NULL);
660 }
661 if (ret == 0) {
662 varName = name + string("_Owner");
663 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
664 (float*)&owner[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
665 NULL);
666 }
667 if (ret == 0) {
668 varName = name + string("_Tag");
669 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
670 (float*)&tag[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
671 NULL);
672 }
673
674 if (reducedElements) {
675 reducedElements->writeToSilo(dbfile, siloPath);
676 }
677
678 // "Elements" is a special case
679 if (name == "Elements") {
680 nodeMesh->writeToSilo(dbfile);
681 }
682
683 return (ret == 0);
684
685 #else // !USE_SILO
686 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 ret.useQuadNodes = false;
699 ret.quadDim = 0;
700
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 ret.useQuadNodes = true;
757 ret.quadDim = 2;
758 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 ret.useQuadNodes = true;
817 ret.quadDim = 3;
818 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 //
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 } // namespace escriptexport
931

  ViewVC Help
Powered by ViewVC 1.1.26