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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2888 - (show annotations)
Fri Jan 29 00:07:00 2010 UTC (9 years, 7 months ago) by caltinay
File size: 27849 byte(s)
Looks like netcdf.hh doesn't exist on Windows so changed to netcdfcpp.h.

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 <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 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 //
571 //
572 //
573 void ElementData::writeConnectivityVTK(ostream& os)
574 {
575 if (numElements > 0) {
576 const IntVec& gNI = nodeMesh->getGlobalNodeIndices();
577 IntVec::const_iterator it;
578 int count = 1;
579 for (it=nodes.begin(); it!=nodes.end(); it++, count++) {
580 os << gNI[*it];
581 if (count % nodesPerElement == 0)
582 os << endl;
583 else
584 os << " ";
585 }
586 }
587 }
588
589 #if USE_SILO
590 //
591 //
592 //
593 inline int toSiloElementType(int type)
594 {
595 switch (type) {
596 case ZONETYPE_BEAM: return DB_ZONETYPE_BEAM;
597 case ZONETYPE_HEX: return DB_ZONETYPE_HEX;
598 case ZONETYPE_POLYGON: return DB_ZONETYPE_POLYGON;
599 case ZONETYPE_QUAD: return DB_ZONETYPE_QUAD;
600 case ZONETYPE_TET: return DB_ZONETYPE_TET;
601 case ZONETYPE_TRIANGLE: return DB_ZONETYPE_TRIANGLE;
602 }
603 return 0;
604 }
605 #endif
606
607 //
608 //
609 //
610 bool ElementData::writeToSilo(DBfile* dbfile, const string& siloPath)
611 {
612 #if USE_SILO
613 if (numElements == 0)
614 return true;
615
616 int ret;
617
618 if (siloPath != "") {
619 ret = DBSetDir(dbfile, siloPath.c_str());
620 if (ret != 0)
621 return false;
622 }
623
624 // write out the full mesh in any case
625 nodeMesh->setSiloPath(siloPath);
626 string siloMeshNameStr = nodeMesh->getFullSiloName();
627 const char* siloMeshName = siloMeshNameStr.c_str();
628 int arraylen = numElements * nodesPerElement;
629 int eltype = toSiloElementType(type);
630
631 string varName = name + string("_zones");
632 ret = DBPutZonelist2(dbfile, varName.c_str(), numElements,
633 nodeMesh->getNumDims(), &nodes[0], arraylen, 0, 0,
634 numGhostElements, &eltype, &nodesPerElement, &numElements, 1, NULL);
635 if (ret == 0) {
636 CoordArray& coordbase = const_cast<CoordArray&>(nodeMesh->getCoords());
637 ret = DBPutUcdmesh(dbfile, siloMeshName,
638 nodeMesh->getNumDims(), NULL, &coordbase[0],
639 nodeMesh->getNumNodes(), numElements, varName.c_str(),
640 /*"facelist"*/NULL, DB_FLOAT, NULL);
641 }
642
643 // Point mesh is useful for debugging
644 if (0) {
645 CoordArray& coordbase = const_cast<CoordArray&>(nodeMesh->getCoords());
646 DBPutPointmesh(dbfile, "/pointmesh", nodeMesh->getNumDims(),
647 &coordbase[0], nodeMesh->getNumNodes(), DB_FLOAT, NULL);
648 }
649
650 if (ret != 0)
651 return false;
652
653 // write out the element-centered variables
654 varName = name + string("_Color");
655 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
656 (float*)&color[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
657 NULL);
658 if (ret == 0) {
659 varName = name + string("_Id");
660 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
661 (float*)&ID[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
662 NULL);
663 }
664 if (ret == 0) {
665 varName = name + string("_Owner");
666 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
667 (float*)&owner[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
668 NULL);
669 }
670 if (ret == 0) {
671 varName = name + string("_Tag");
672 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
673 (float*)&tag[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
674 NULL);
675 }
676
677 if (reducedElements) {
678 reducedElements->writeToSilo(dbfile, siloPath);
679 }
680
681 // "Elements" is a special case
682 if (name == "Elements") {
683 nodeMesh->writeToSilo(dbfile);
684 }
685
686 return (ret == 0);
687
688 #else // !USE_SILO
689 return false;
690 #endif
691 }
692
693 //
694 //
695 //
696 FinleyElementInfo ElementData::getFinleyTypeInfo(ElementTypeId typeId)
697 {
698 FinleyElementInfo ret;
699 ret.multiCellIndices = NULL;
700 ret.elementFactor = 1;
701 ret.useQuadNodes = false;
702 ret.quadDim = 0;
703
704 switch (typeId) {
705 case Point1_Contact://untested
706 case Line2Face_Contact://untested
707 case Line3Face_Contact://untested
708 case Line2Face://untested
709 case Line3Face://untested
710 case Point1://untested
711 cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
712 ret.elementSize = 1;
713 ret.elementType = ZONETYPE_POLYGON;
714 break;
715
716 case Tri3Face_Contact://untested
717 case Tri3Face://untested
718 cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
719 case Line2_Contact:
720 case Rec4Face_Contact:
721 case Rec4Face:
722 case Line2:
723 ret.elementSize = ret.reducedElementSize = 2;
724 ret.elementType = ret.reducedElementType = ZONETYPE_BEAM;
725 break;
726
727 case Line3:
728 ret.multiCellIndices = line3indices;
729 ret.elementFactor = 2;
730 // fall through
731 case Line3_Contact:
732 case Tri6Face_Contact://untested
733 case Rec8Face_Contact:
734 case Tri6Face://untested
735 case Rec8Face:
736 //VTK_QUADRATIC_EDGE
737 ret.elementSize = ret.reducedElementSize = 2;
738 ret.elementType = ret.reducedElementType = ZONETYPE_BEAM;
739 break;
740
741 case Tet4Face_Contact://untested
742 case Tet4Face://untested
743 cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
744 case Tri3_Contact:
745 case Tri3:
746 ret.elementSize = ret.reducedElementSize = 3;
747 ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
748 break;
749
750 case Rec4_Contact:
751 case Hex8Face_Contact:
752 case Hex8Face:
753 case Rec4:
754 ret.elementSize = ret.reducedElementSize = 4;
755 ret.elementType = ret.reducedElementType = ZONETYPE_QUAD;
756 break;
757
758 case Rec9:
759 ret.useQuadNodes = true;
760 ret.quadDim = 2;
761 ret.multiCellIndices = rec9indices;
762 ret.elementFactor = 4;
763 // fall through
764 case Rec9_Contact:
765 ret.elementSize = ret.reducedElementSize = 4;
766 ret.elementType = ret.reducedElementType = ZONETYPE_QUAD;
767 break;
768
769 case Tet4:
770 ret.elementSize = ret.reducedElementSize = 4;
771 ret.elementType = ret.reducedElementType = ZONETYPE_TET;
772 break;
773
774 case Tri6:
775 ret.multiCellIndices = tri6indices;
776 ret.elementFactor = 4;
777 // fall through
778 case Tri6_Contact:
779 case Tet10Face_Contact://untested
780 case Tet10Face://untested
781 //VTK_QUADRATIC_TRIANGLE
782 ret.elementSize = ret.reducedElementSize = 3;
783 ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
784 break;
785
786 case Rec8:
787 ret.multiCellIndices = rec8indices;
788 ret.elementFactor = 6;
789 // fall through
790 case Hex20Face:
791 case Rec8_Contact:
792 case Hex20Face_Contact:
793 //VTK_QUADRATIC_QUAD
794 ret.elementSize = 3;
795 ret.elementType = ZONETYPE_TRIANGLE;
796 ret.reducedElementSize = 4;
797 ret.reducedElementType = ZONETYPE_QUAD;
798 break;
799
800 case Tet10:
801 //VTK_QUADRATIC_TETRA
802 ret.multiCellIndices = tet10indices;
803 ret.elementFactor = 8;
804 ret.elementSize = ret.reducedElementSize = 4;
805 ret.elementType = ret.reducedElementType = ZONETYPE_TET;
806 break;
807
808 case Hex20:
809 //VTK_QUADRATIC_HEXAHEDRON
810 ret.multiCellIndices = hex20indices;
811 ret.elementFactor = 36;
812 ret.elementSize = 3;
813 ret.elementType = ZONETYPE_TRIANGLE;
814 ret.reducedElementSize = 8;
815 ret.reducedElementType = ZONETYPE_HEX;
816 break;
817
818 case Hex27:
819 ret.useQuadNodes = true;
820 ret.quadDim = 3;
821 ret.multiCellIndices = hex27indices;
822 ret.elementFactor = 8;
823 // fall through
824 case Hex8:
825 ret.elementSize = ret.reducedElementSize = 8;
826 ret.elementType = ret.reducedElementType = ZONETYPE_HEX;
827 break;
828
829 default:
830 cerr << "WARNING: Unknown Finley Type " << typeId << endl;
831 break;
832 }
833 return ret;
834 }
835
836 //
837 //
838 //
839 inline bool inside1D(float x, float c, float r)
840 {
841 return (ABS(x-c) <= r);
842 }
843
844 //
845 //
846 //
847 inline bool inside2D(float x, float y, float cx, float cy, float r)
848 {
849 return (inside1D(x, cx, r) && inside1D(y, cy, r));
850 }
851
852 //
853 //
854 //
855 inline bool inside3D(float x, float y, float z,
856 float cx, float cy, float cz, float r)
857 {
858 return (inside2D(x, y, cx, cy, r) && inside1D(z, cz, r));
859 }
860
861 //
862 //
863 //
864 QuadMaskInfo ElementData::buildQuadMask(const CoordArray& qnodes, int numQNodes)
865 {
866 QuadMaskInfo qmi;
867
868 if (numQNodes == 0)
869 return qmi;
870
871 if (finleyTypeId == Line3Macro) {
872 for (int i=0; i<elementFactor; i++) {
873 const float bounds[] = { 0.25, 0.75 };
874 IntVec m(numQNodes, 0);
875 int hits = 0;
876 for (size_t j=0; j<numQNodes; j++) {
877 if (inside1D(qnodes[0][j], bounds[i], .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 == Rec9) || (finleyTypeId == Rec9Macro)) {
889 for (int i=0; i<elementFactor; i++) {
890 const float bounds[][2] = { { 0.25, 0.25 }, { 0.75, 0.25 },
891 { 0.25, 0.75 }, { 0.75, 0.75 } };
892 IntVec m(numQNodes, 0);
893 int hits = 0;
894 for (size_t j=0; j<numQNodes; j++) {
895 if (inside2D(qnodes[0][j], qnodes[1][j],
896 bounds[i][0], bounds[i][1], 0.25)) {
897 m[j] = 1;
898 hits++;
899 }
900 }
901 qmi.mask.push_back(m);
902 if (hits == 0)
903 qmi.factor.push_back(1);
904 else
905 qmi.factor.push_back(hits);
906 }
907 } else if ((finleyTypeId == Hex27) || (finleyTypeId == Hex27Macro) ){
908 for (int i=0; i<elementFactor; i++) {
909 const float bounds[][3] = {
910 { 0.25, 0.25, 0.25 }, { 0.75, 0.25, 0.25 },
911 { 0.25, 0.75, 0.25 }, { 0.75, 0.75, 0.25 },
912 { 0.25, 0.25, 0.75 }, { 0.75, 0.25, 0.75 },
913 { 0.25, 0.75, 0.75 }, { 0.75, 0.75, 0.75 } };
914 IntVec m(numQNodes, 0);
915 int hits = 0;
916 for (size_t j=0; j<numQNodes; j++) {
917 if (inside3D(qnodes[0][j], qnodes[1][j], qnodes[2][j],
918 bounds[i][0], bounds[i][1], bounds[i][2], 0.25)) {
919 m[j] = 1;
920 hits++;
921 }
922 }
923 qmi.mask.push_back(m);
924 if (hits == 0)
925 qmi.factor.push_back(1);
926 else
927 qmi.factor.push_back(hits);
928 }
929 }
930 return qmi;
931 }
932
933 } // namespace escriptexport
934

  ViewVC Help
Powered by ViewVC 1.1.26