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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2940 - (show annotations)
Fri Feb 19 00:38:45 2010 UTC (9 years, 6 months ago) by caltinay
File size: 34663 byte(s)
Minor potential optimization: prefer empty() over size()==0.

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 5, 4, 2,
43 3, 1, 4,
44 4, 5, 3
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, 4, 0, 7,
62 6, 5, 4, 8,
63 5, 1, 4, 8,
64 9, 8, 7, 3,
65 2, 5, 6, 9,
66 8, 9, 5, 6,
67 6, 7, 9, 8,
68 6, 4, 7, 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.empty())
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 case Line3Macro:
737 ret.multiCellIndices = line3indices;
738 ret.elementFactor = 2;
739 // fall through
740 case Line3_Contact:
741 case Tri6Face_Contact://untested
742 case Rec8Face_Contact:
743 case Tri6Face://untested
744 case Rec8Face:
745 //VTK_QUADRATIC_EDGE
746 ret.elementSize = ret.reducedElementSize = 2;
747 ret.elementType = ret.reducedElementType = ZONETYPE_BEAM;
748 break;
749
750 case Tet4Face_Contact://untested
751 case Tet4Face://untested
752 cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
753 case Tri3_Contact:
754 case Tri3:
755 ret.elementSize = ret.reducedElementSize = 3;
756 ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
757 break;
758
759 case Rec4_Contact:
760 case Hex8Face_Contact:
761 case Hex8Face:
762 case Rec4:
763 ret.elementSize = ret.reducedElementSize = 4;
764 ret.elementType = ret.reducedElementType = ZONETYPE_QUAD;
765 break;
766
767 case Rec9:
768 case Rec9Macro:
769 ret.useQuadNodes = true;
770 ret.quadDim = 2;
771 ret.multiCellIndices = rec9indices;
772 ret.elementFactor = 4;
773 // fall through
774 case Rec9_Contact:
775 ret.elementSize = ret.reducedElementSize = 4;
776 ret.elementType = ret.reducedElementType = ZONETYPE_QUAD;
777 break;
778
779 case Tet4:
780 ret.elementSize = ret.reducedElementSize = 4;
781 ret.elementType = ret.reducedElementType = ZONETYPE_TET;
782 break;
783
784 case Tri6:
785 case Tri6Macro:
786 ret.useQuadNodes = true;
787 ret.quadDim = 2;
788 ret.multiCellIndices = tri6indices;
789 ret.elementFactor = 4;
790 // fall through
791 case Tri6_Contact:
792 case Tet10Face_Contact://untested
793 case Tet10Face://untested
794 //VTK_QUADRATIC_TRIANGLE
795 ret.elementSize = ret.reducedElementSize = 3;
796 ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
797 break;
798
799 case Rec8:
800 ret.multiCellIndices = rec8indices;
801 ret.elementFactor = 6;
802 // fall through
803 case Hex20Face:
804 case Rec8_Contact:
805 case Hex20Face_Contact:
806 //VTK_QUADRATIC_QUAD
807 ret.elementSize = 3;
808 ret.elementType = ZONETYPE_TRIANGLE;
809 ret.reducedElementSize = 4;
810 ret.reducedElementType = ZONETYPE_QUAD;
811 break;
812
813 case Tet10:
814 case Tet10Macro:
815 //VTK_QUADRATIC_TETRA
816 ret.useQuadNodes = true;
817 ret.quadDim = 3;
818 ret.multiCellIndices = tet10indices;
819 ret.elementFactor = 8;
820 ret.elementSize = ret.reducedElementSize = 4;
821 ret.elementType = ret.reducedElementType = ZONETYPE_TET;
822 break;
823
824 case Hex20:
825 //VTK_QUADRATIC_HEXAHEDRON
826 ret.multiCellIndices = hex20indices;
827 ret.elementFactor = 36;
828 ret.elementSize = 3;
829 ret.elementType = ZONETYPE_TRIANGLE;
830 ret.reducedElementSize = 8;
831 ret.reducedElementType = ZONETYPE_HEX;
832 break;
833
834 case Hex27:
835 case Hex27Macro:
836 ret.useQuadNodes = true;
837 ret.quadDim = 3;
838 ret.multiCellIndices = hex27indices;
839 ret.elementFactor = 8;
840 // fall through
841 case Hex8:
842 ret.elementSize = ret.reducedElementSize = 8;
843 ret.elementType = ret.reducedElementType = ZONETYPE_HEX;
844 break;
845
846 default:
847 cerr << "WARNING: Unknown Finley Type " << typeId << endl;
848 break;
849 }
850 return ret;
851 }
852
853 /////////////////////////////////
854 // Helpers for buildQuadMask() //
855 /////////////////////////////////
856
857 // returns true if |x-c| <= r, false otherwise
858 inline bool inside1D(float x, float c, float r)
859 {
860 return (ABS(x-c) <= r);
861 }
862
863 // returns true if |x-cx| <= r and |y-cy| <= r, false otherwise
864 inline bool inside2D(float x, float y, float cx, float cy, float r)
865 {
866 return (inside1D(x, cx, r) && inside1D(y, cy, r));
867 }
868
869 // returns true if |x-cx| <= r and |y-cy| <= r and |z-cz| <= r, false otherwise
870 inline bool inside3D(float x, float y, float z,
871 float cx, float cy, float cz, float r)
872 {
873 return (inside2D(x, y, cx, cy, r) && inside1D(z, cz, r));
874 }
875
876 // returns true if d1 and d2 have the same sign or at least one of them is
877 // close to 0, false otherwise
878 inline bool sameSide(float d1, float d2)
879 {
880 const float TOL = 1.e-8f;
881 return (ABS(d1) < TOL || ABS(d2) < TOL || d1*d2>=0.);
882 }
883
884 // computes the determinant of the 4x4 matrix given by its elements m_ij
885 static float det4x4(float m_00, float m_01, float m_02, float m_03,
886 float m_10, float m_11, float m_12, float m_13,
887 float m_20, float m_21, float m_22, float m_23,
888 float m_30, float m_31, float m_32, float m_33)
889 {
890 float det1 = m_12 * m_23 - m_22 * m_13;
891 float det2 = m_11 * m_23 - m_21 * m_13;
892 float det3 = m_11 * m_22 - m_21 * m_12;
893 float det4 = m_10 * m_23 - m_20 * m_13;
894 float det5 = m_10 * m_22 - m_20 * m_12;
895 float det6 = m_10 * m_21 - m_20 * m_11;
896 return -m_30 * (m_01 * det1 - m_02 * det2 + m_03 * det3) +
897 m_31 * (m_00 * det1 - m_02 * det4 + m_03 * det5) -
898 m_32 * (m_00 * det2 - m_01 * det4 + m_03 * det6) +
899 m_33 * (m_00 * det3 - m_01 * det5 + m_02 * det6);
900 }
901
902 // returns true if point (x,y,z) is in or on the tetrahedron given by its
903 // corner points p0, p1, p2 and p3, false otherwise.
904 static bool pointInTet(float x, float y, float z,
905 const float* p0, const float* p1,
906 const float* p2, const float* p3)
907 {
908 float d0 = det4x4(
909 p0[0], p0[1], p0[2], 1.f,
910 p1[0], p1[1], p1[2], 1.f,
911 p2[0], p2[1], p2[2], 1.f,
912 p3[0], p3[1], p3[2], 1.f);
913 float d1 = det4x4(
914 x, y, z, 1.f,
915 p1[0], p1[1], p1[2], 1.f,
916 p2[0], p2[1], p2[2], 1.f,
917 p3[0], p3[1], p3[2], 1.f);
918 float d2 = det4x4(
919 p0[0], p0[1], p0[2], 1.f,
920 x, y, z, 1.f,
921 p2[0], p2[1], p2[2], 1.f,
922 p3[0], p3[1], p3[2], 1.f);
923 float d3 = det4x4(
924 p0[0], p0[1], p0[2], 1.f,
925 p1[0], p1[1], p1[2], 1.f,
926 x, y, z, 1.f,
927 p3[0], p3[1], p3[2], 1.f);
928 float d4 = det4x4(
929 p0[0], p0[1], p0[2], 1.f,
930 p1[0], p1[1], p1[2], 1.f,
931 p2[0], p2[1], p2[2], 1.f,
932 x, y, z, 1.f);
933
934 return (sameSide(d0,d1) && sameSide(d1,d2) &&
935 sameSide(d2,d3) && sameSide(d3,d4));
936 }
937
938 // returns true if point (x,y) is in or on the triangle given by its corner
939 // points p0, p1 and p2, false otherwise
940 static bool pointInTri(float x, float y,
941 const float* p0, const float* p1, const float* p2)
942 {
943 const float TOL = 1.e-8f;
944 float v0[2] = { p2[0]-p0[0], p2[1]-p0[1] };
945 float v1[2] = { p1[0]-p0[0], p1[1]-p0[1] };
946 float v2[2] = { x - p0[0], y - p0[1] };
947
948 float dot00 = v0[0]*v0[0]+v0[1]*v0[1];
949 float dot01 = v0[0]*v1[0]+v0[1]*v1[1];
950 float dot02 = v0[0]*v2[0]+v0[1]*v2[1];
951 float dot11 = v1[0]*v1[0]+v1[1]*v1[1];
952 float dot12 = v1[0]*v2[0]+v1[1]*v2[1];
953 float invDenom = dot00*dot11 - dot01*dot01;
954 if (ABS(invDenom) < TOL) invDenom = TOL;
955 invDenom = 1./invDenom;
956 float u = (dot11*dot02 - dot01*dot12) * invDenom;
957 float v = (dot00*dot12 - dot01*dot02) * invDenom;
958 return (u>=0.f) && (v>=0.f) && (u+v<=1.f);
959 }
960
961 //
962 //
963 //
964 QuadMaskInfo ElementData::buildQuadMask(const CoordArray& qnodes, int numQNodes)
965 {
966 QuadMaskInfo qmi;
967 if (numQNodes == 0)
968 return qmi;
969
970 if (finleyTypeId == Line3Macro) {
971 for (int i=0; i<elementFactor; i++) {
972 const float bounds[] = { 0.25, 0.75 };
973 IntVec m(numQNodes, 0);
974 int hits = 0;
975 for (size_t j=0; j<numQNodes; j++) {
976 if (inside1D(qnodes[0][j], bounds[i], .25)) {
977 m[j] = 1;
978 hits++;
979 }
980 }
981 qmi.mask.push_back(m);
982 if (hits == 0)
983 qmi.factor.push_back(1);
984 else
985 qmi.factor.push_back(hits);
986 }
987 } else if ((finleyTypeId == Tri6) || (finleyTypeId == Tri6Macro)) {
988 for (int i=0; i<elementFactor; i++) {
989 const float bounds[][2] = { { 0., 0. }, { 1., 0. },
990 { 0., 1. }, { .5, 0. },
991 { .5, .5 }, { 0., .5 } };
992 const size_t* nodeIdx = &tri6indices[i*nodesPerElement];
993 IntVec m(numQNodes, 0);
994 int hits = 0;
995 for (size_t j=0; j<numQNodes; j++) {
996 // check if point j is in triangle i
997 if (pointInTri(qnodes[0][j], qnodes[1][j],
998 bounds[nodeIdx[0]], bounds[nodeIdx[1]],
999 bounds[nodeIdx[2]])) {
1000 m[j] = 1;
1001 hits++;
1002 }
1003 }
1004 if (hits == 0) {
1005 // if an element does not contain any quadrature points we
1006 // simply average over all data points within that element
1007 m = IntVec(numQNodes, 1);
1008 qmi.factor.push_back(numQNodes);
1009 qmi.factor.push_back(1);
1010 } else {
1011 qmi.factor.push_back(hits);
1012 }
1013 qmi.mask.push_back(m);
1014 }
1015 } else if ((finleyTypeId == Tet10) || (finleyTypeId == Tet10Macro)) {
1016 for (int i=0; i<elementFactor; i++) {
1017 const float bounds[][3] = {
1018 { 0., 0., 0. },
1019 { 1., 0., 0. },
1020 { 0., 0., 1. },
1021 { 0., 1., 0. },
1022 { .5, 0., 0. },
1023 { .5, 0., .5 },
1024 { 0., 0., .5 },
1025 { 0., .5, 0. },
1026 { .5, .5, 0. },
1027 { 0., .5, .5 }
1028 };
1029 // need to reorder the elements
1030 const size_t elNumIdx[] = { 0,1,2,4,3,5,6,7 };
1031 const size_t* nodeIdx = &tet10indices[elNumIdx[i]*nodesPerElement];
1032 IntVec m(numQNodes, 0);
1033 int hits = 0;
1034 for (size_t j=0; j<numQNodes; j++) {
1035 // check if point j is in tetrahedron i
1036 if (pointInTet(qnodes[0][j], qnodes[1][j], qnodes[2][j],
1037 bounds[nodeIdx[0]], bounds[nodeIdx[1]],
1038 bounds[nodeIdx[2]], bounds[nodeIdx[3]])) {
1039 m[j] = 1;
1040 hits++;
1041 }
1042 }
1043 if (hits == 0) {
1044 // if an element does not contain any quadrature points we
1045 // simply average over all data points within that element
1046 m = IntVec(numQNodes, 1);
1047 qmi.factor.push_back(numQNodes);
1048 } else {
1049 qmi.factor.push_back(hits);
1050 }
1051 qmi.mask.push_back(m);
1052 }
1053 } else if ((finleyTypeId == Rec9) || (finleyTypeId == Rec9Macro)) {
1054 for (int i=0; i<elementFactor; i++) {
1055 const float bounds[][2] = { { 0.25, 0.25 }, { 0.75, 0.25 },
1056 { 0.25, 0.75 }, { 0.75, 0.75 } };
1057 IntVec m(numQNodes, 0);
1058 int hits = 0;
1059 for (size_t j=0; j<numQNodes; j++) {
1060 if (inside2D(qnodes[0][j], qnodes[1][j],
1061 bounds[i][0], bounds[i][1], 0.25)) {
1062 m[j] = 1;
1063 hits++;
1064 }
1065 }
1066 qmi.mask.push_back(m);
1067 if (hits == 0)
1068 qmi.factor.push_back(1);
1069 else
1070 qmi.factor.push_back(hits);
1071 }
1072 } else if ((finleyTypeId == Hex27) || (finleyTypeId == Hex27Macro) ){
1073 for (int i=0; i<elementFactor; i++) {
1074 const float bounds[][3] = {
1075 { 0.25, 0.25, 0.25 }, { 0.75, 0.25, 0.25 },
1076 { 0.25, 0.75, 0.25 }, { 0.75, 0.75, 0.25 },
1077 { 0.25, 0.25, 0.75 }, { 0.75, 0.25, 0.75 },
1078 { 0.25, 0.75, 0.75 }, { 0.75, 0.75, 0.75 } };
1079 IntVec m(numQNodes, 0);
1080 int hits = 0;
1081 for (size_t j=0; j<numQNodes; j++) {
1082 if (inside3D(qnodes[0][j], qnodes[1][j], qnodes[2][j],
1083 bounds[i][0], bounds[i][1], bounds[i][2], 0.25)) {
1084 m[j] = 1;
1085 hits++;
1086 }
1087 }
1088 qmi.mask.push_back(m);
1089 if (hits == 0)
1090 qmi.factor.push_back(1);
1091 else
1092 qmi.factor.push_back(hits);
1093 }
1094 }
1095 return qmi;
1096 }
1097
1098 } // namespace escriptexport
1099

  ViewVC Help
Powered by ViewVC 1.1.26