/[escript]/trunk/weipa/src/FinleyElements.cpp
ViewVC logotype

Contents of /trunk/weipa/src/FinleyElements.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3398 - (show annotations)
Mon Dec 6 00:02:33 2010 UTC (10 years, 2 months ago) by caltinay
File size: 39093 byte(s)
Fixed the quad mask in a very specific case and changed saveVTK to write a
specific value for unused nodes instead of the first one which breaks unit
tests depending on the number of ranks used.
Still need to rebaseline some VTK files so the tests pass.

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 <weipa/FinleyElements.h>
15 #include <weipa/NodeData.h>
16
17 #ifndef VISIT_PLUGIN
18 #include <dudley/CppAdapter/MeshAdapter.h>
19 #include <finley/CppAdapter/MeshAdapter.h>
20 #else
21 #define ABS(X) ((X)>0?(X):-(X))
22 #endif
23
24 #include <iostream>
25
26 #if USE_NETCDF
27 #include <netcdfcpp.h>
28 #endif
29
30 #if USE_SILO
31 #include <silo.h>
32 #endif
33
34
35 using namespace std;
36
37 // The following arrays contain indices to convert unsupported element
38 // types into supported ones (lines, triangles, quads, hexahedrons)
39
40 static const size_t line3indices[2*2] = {
41 0, 2,
42 2, 1
43 };
44 static const size_t tri6indices[4*3] = {
45 0, 3, 5,
46 5, 4, 2,
47 3, 1, 4,
48 4, 5, 3
49 };
50 static const size_t rec8indices[6*3] = {
51 0, 4, 7,
52 4, 1, 5,
53 5, 2, 6,
54 6, 3, 7,
55 7, 5, 6,
56 7, 4, 5
57 };
58 static const size_t rec9indices[4*4] = {
59 0, 4, 8, 7,
60 4, 1, 5, 8,
61 7, 8, 6, 3,
62 8, 5, 2, 6
63 };
64 static const size_t tet10indices[8*4] = {
65 6, 4, 0, 7,
66 6, 5, 4, 8,
67 5, 1, 4, 8,
68 9, 8, 7, 3,
69 2, 5, 6, 9,
70 8, 9, 5, 6,
71 6, 7, 9, 8,
72 6, 4, 7, 8
73 };
74 static const size_t hex20indices[36*3] = {
75 0, 8, 12, 8, 1, 13, 13, 5, 16,
76 16, 4, 12, 8, 13, 16, 8, 16, 12,
77 1, 9, 13, 9, 2, 14, 14, 6, 17,
78 17, 5, 13, 9, 14, 17, 9, 17, 13,
79 2, 10, 14, 10, 3, 15, 15, 7, 18,
80 18, 14, 6, 10, 15, 18, 10, 18, 14,
81 3, 11, 15, 11, 0, 12, 12, 4, 19,
82 19, 7, 15, 11, 12, 19, 11, 19, 15,
83 4, 16, 19, 16, 5, 17, 17, 6, 18,
84 18, 7, 19, 16, 17, 18, 16, 18, 19,
85 3, 10, 11, 10, 2, 9, 9, 1, 8,
86 8, 0, 11, 10, 9, 8, 10, 8, 11
87 };
88 static const size_t hex27indices[8*8] = {
89 0, 8, 20, 11, 12, 21, 26, 24,
90 8, 1, 9, 20, 21, 13, 22, 26,
91 11, 20, 10, 3, 24, 26, 23, 15,
92 20, 9, 2, 10, 26, 22, 14, 23,
93 12, 21, 26, 24, 4, 16, 25, 19,
94 21, 13, 22, 26, 16, 5, 17, 25,
95 24, 26, 23, 15, 19, 25, 18, 7,
96 26, 22, 14, 23, 25, 17, 6, 18
97 };
98
99 namespace weipa {
100
101 //
102 // Constructor
103 //
104 FinleyElements::FinleyElements(const string& elementName, FinleyNodes_ptr nodeData)
105 : originalMesh(nodeData), name(elementName), numElements(0),
106 numGhostElements(0), nodesPerElement(0),
107 type(ZONETYPE_UNKNOWN), finleyTypeId(Finley_NoRef), elementFactor(1)
108 {
109 nodeMesh.reset(new FinleyNodes(name));
110 }
111
112 //
113 // Copy constructor
114 //
115 FinleyElements::FinleyElements(const FinleyElements& e)
116 {
117 name = e.name;
118 numElements = e.numElements;
119 numGhostElements = e.numGhostElements;
120 type = e.type;
121 finleyTypeId = e.finleyTypeId;
122 nodesPerElement = e.nodesPerElement;
123 elementFactor = e.elementFactor;
124 originalMesh = e.originalMesh;
125 if (e.nodeMesh)
126 nodeMesh.reset(new FinleyNodes(*e.nodeMesh));
127 else
128 nodeMesh.reset(new FinleyNodes(name));
129
130 nodes = e.nodes;
131 ID = e.ID;
132 color = e.color;
133 tag = e.tag;
134 owner = e.owner;
135
136 if (e.reducedElements)
137 reducedElements = FinleyElements_ptr(
138 new FinleyElements(*e.reducedElements));
139 }
140
141 //
142 //
143 //
144 bool FinleyElements::initFromDudley(const Dudley_ElementFile* dudleyFile)
145 {
146 #ifndef VISIT_PLUGIN
147 numElements = dudleyFile->numElements;
148
149 if (numElements > 0) {
150 nodesPerElement = dudleyFile->numNodes;
151
152 int* iPtr;
153
154 iPtr = dudleyFile->Nodes;
155 nodes.clear();
156 nodes.insert(nodes.end(), numElements*nodesPerElement, 0);
157 copy(iPtr, iPtr+numElements*nodesPerElement, nodes.begin());
158
159 iPtr = dudleyFile->Color;
160 color.clear();
161 color.insert(color.end(), numElements, 0);
162 copy(iPtr, iPtr+numElements, color.begin());
163
164 iPtr = dudleyFile->Id;
165 ID.clear();
166 ID.insert(ID.end(), numElements, 0);
167 copy(iPtr, iPtr+numElements, ID.begin());
168
169 iPtr = dudleyFile->Owner;
170 owner.clear();
171 owner.insert(owner.end(), numElements, 0);
172 copy(iPtr, iPtr+numElements, owner.begin());
173
174 iPtr = dudleyFile->Tag;
175 tag.clear();
176 tag.insert(tag.end(), numElements, 0);
177 copy(iPtr, iPtr+numElements, tag.begin());
178
179 FinleyElementInfo f = getDudleyTypeInfo(dudleyFile->etype);
180 type = f.elementType;
181 elementFactor = f.elementFactor;
182 if (elementFactor > 1 || f.reducedElementSize != nodesPerElement)
183 buildReducedElements(f);
184
185 buildMeshes();
186 }
187 return true;
188
189 #else // VISIT_PLUGIN
190 return false;
191 #endif
192 }
193
194 //
195 //
196 //
197 bool FinleyElements::initFromFinley(const Finley_ElementFile* finleyFile)
198 {
199 #ifndef VISIT_PLUGIN
200 numElements = finleyFile->numElements;
201
202 if (numElements > 0) {
203 nodesPerElement = finleyFile->numNodes;
204
205 int* iPtr;
206
207 iPtr = finleyFile->Nodes;
208 nodes.clear();
209 nodes.insert(nodes.end(), numElements*nodesPerElement, 0);
210 copy(iPtr, iPtr+numElements*nodesPerElement, nodes.begin());
211
212 iPtr = finleyFile->Color;
213 color.clear();
214 color.insert(color.end(), numElements, 0);
215 copy(iPtr, iPtr+numElements, color.begin());
216
217 iPtr = finleyFile->Id;
218 ID.clear();
219 ID.insert(ID.end(), numElements, 0);
220 copy(iPtr, iPtr+numElements, ID.begin());
221
222 iPtr = finleyFile->Owner;
223 owner.clear();
224 owner.insert(owner.end(), numElements, 0);
225 copy(iPtr, iPtr+numElements, owner.begin());
226
227 iPtr = finleyFile->Tag;
228 tag.clear();
229 tag.insert(tag.end(), numElements, 0);
230 copy(iPtr, iPtr+numElements, tag.begin());
231
232 finleyTypeId = finleyFile->referenceElementSet->referenceElement
233 ->Type->TypeId;
234 FinleyElementInfo f = getFinleyTypeInfo(finleyTypeId);
235 type = f.elementType;
236 elementFactor = f.elementFactor;
237 if (elementFactor > 1 || f.reducedElementSize != nodesPerElement)
238 buildReducedElements(f);
239
240 if (f.useQuadNodes) {
241 CoordArray quadNodes;
242 int numQuadNodes;
243 Finley_ShapeFunction* sf = finleyFile->referenceElementSet
244 ->referenceElement->Parametrization;
245 numQuadNodes = sf->numQuadNodes;
246 for (int i=0; i<f.quadDim; i++) {
247 double* srcPtr = sf->QuadNodes+i;
248 float* c = new float[numQuadNodes];
249 quadNodes.push_back(c);
250 for (int j=0; j<numQuadNodes; j++, srcPtr+=f.quadDim) {
251 *c++ = (float) *srcPtr;
252 }
253 }
254 quadMask = buildQuadMask(quadNodes, numQuadNodes);
255 for (int i=0; i<f.quadDim; i++)
256 delete[] quadNodes[i];
257 quadNodes.clear();
258
259
260 // now the reduced quadrature
261 sf = finleyFile->referenceElementSet
262 ->referenceElementReducedQuadrature->Parametrization;
263 numQuadNodes = sf->numQuadNodes;
264 for (int i=0; i<f.quadDim; i++) {
265 double* srcPtr = sf->QuadNodes+i;
266 float* c = new float[numQuadNodes];
267 quadNodes.push_back(c);
268 for (int j=0; j<numQuadNodes; j++, srcPtr+=f.quadDim) {
269 *c++ = (float) *srcPtr;
270 }
271 }
272 reducedQuadMask = buildQuadMask(quadNodes, numQuadNodes);
273 for (int i=0; i<f.quadDim; i++)
274 delete[] quadNodes[i];
275 quadNodes.clear();
276 }
277
278 buildMeshes();
279 }
280 return true;
281
282 #else // VISIT_PLUGIN
283 return false;
284 #endif
285 }
286
287 //
288 // Reads element data from given NetCDF file
289 //
290 bool FinleyElements::readFromNc(NcFile* ncfile)
291 {
292 #if USE_NETCDF
293 string num_str("num_");
294 num_str += name;
295
296 NcAtt* att = ncfile->get_att(num_str.c_str());
297 numElements = att->as_int(0);
298
299 // Only attempt to read further if there are any elements.
300 // Having no elements is not an error.
301 if (numElements > 0) {
302 att = ncfile->get_att((num_str + string("_numNodes")).c_str());
303 nodesPerElement = att->as_int(0);
304
305 nodes.insert(nodes.end(), numElements*nodesPerElement, 0);
306 NcVar* var = ncfile->get_var((name + string("_Nodes")).c_str());
307 var->get(&nodes[0], numElements, nodesPerElement);
308
309 color.insert(color.end(), numElements, 0);
310 var = ncfile->get_var((name + string("_Color")).c_str());
311 var->get(&color[0], numElements);
312
313 ID.insert(ID.end(), numElements, 0);
314 var = ncfile->get_var((name + string("_Id")).c_str());
315 var->get(&ID[0], numElements);
316
317 owner.insert(owner.end(), numElements, 0);
318 var = ncfile->get_var((name + string("_Owner")).c_str());
319 var->get(&owner[0], numElements);
320
321 tag.insert(tag.end(), numElements, 0);
322 var = ncfile->get_var((name + string("_Tag")).c_str());
323 var->get(&tag[0], numElements);
324
325 att = ncfile->get_att((name + string("_TypeId")).c_str());
326 finleyTypeId = (Finley_ElementTypeId)att->as_int(0);
327 FinleyElementInfo f = getFinleyTypeInfo(finleyTypeId);
328 type = f.elementType;
329 elementFactor = f.elementFactor;
330 if (f.elementFactor > 1 || f.reducedElementSize != nodesPerElement)
331 buildReducedElements(f);
332
333 // if we don't link with finley we can't get the quadrature nodes
334 // and hence cannot interpolate data properly
335 #ifndef VISIT_PLUGIN
336 if (f.useQuadNodes) {
337 att = ncfile->get_att("order");
338 int order = att->as_int(0);
339 att = ncfile->get_att("reduced_order");
340 int reduced_order = att->as_int(0);
341 finley::ReferenceElementSetWrapper wrapper(finleyTypeId, order,
342 reduced_order);
343 Finley_ReferenceElementSet* refElements = wrapper.getElementSet();
344
345 CoordArray quadNodes;
346 int numQuadNodes;
347 Finley_ShapeFunction* sf = refElements->referenceElement
348 ->Parametrization;
349 numQuadNodes = sf->numQuadNodes;
350 for (int i=0; i<f.quadDim; i++) {
351 double* srcPtr = sf->QuadNodes+i;
352 float* c = new float[numQuadNodes];
353 quadNodes.push_back(c);
354 for (int j=0; j<numQuadNodes; j++, srcPtr+=f.quadDim) {
355 *c++ = (float) *srcPtr;
356 }
357 }
358 quadMask = buildQuadMask(quadNodes, numQuadNodes);
359 for (int i=0; i<f.quadDim; i++)
360 delete[] quadNodes[i];
361 quadNodes.clear();
362
363 // now the reduced quadrature
364 sf = refElements->referenceElementReducedQuadrature->Parametrization;
365 numQuadNodes = sf->numQuadNodes;
366 for (int i=0; i<f.quadDim; i++) {
367 double* srcPtr = sf->QuadNodes+i;
368 float* c = new float[numQuadNodes];
369 quadNodes.push_back(c);
370 for (int j=0; j<numQuadNodes; j++, srcPtr+=f.quadDim) {
371 *c++ = (float) *srcPtr;
372 }
373 }
374 reducedQuadMask = buildQuadMask(quadNodes, numQuadNodes);
375 for (int i=0; i<f.quadDim; i++)
376 delete[] quadNodes[i];
377 quadNodes.clear();
378 }
379 #endif // VISIT_PLUGIN
380
381 buildMeshes();
382 }
383
384 return true;
385 #else // !USE_NETCDF
386 return false;
387 #endif
388 }
389
390 //
391 //
392 //
393 StringVec FinleyElements::getMeshNames() const
394 {
395 StringVec res;
396 if (nodeMesh)
397 res.push_back(nodeMesh->getName());
398 if (reducedElements) {
399 StringVec rNames = reducedElements->getMeshNames();
400 if (!rNames.empty())
401 res.insert(res.end(), rNames.begin(), rNames.end());
402 }
403 return res;
404 }
405
406 //
407 //
408 //
409 StringVec FinleyElements::getVarNames() const
410 {
411 StringVec res;
412 res.push_back(name + string("_Color"));
413 res.push_back(name + string("_Id"));
414 res.push_back(name + string("_Owner"));
415 res.push_back(name + string("_Tag"));
416 return res;
417 }
418
419 //
420 //
421 //
422 const IntVec& FinleyElements::getVarDataByName(const string varName) const
423 {
424 if (varName == name+string("_Color"))
425 return color;
426 else if (varName == name+string("_Id"))
427 return ID;
428 else if (varName == name+string("_Owner")) {
429 return owner;
430 } else if (varName == name+string("_Tag"))
431 return tag;
432 else if (reducedElements)
433 return reducedElements->getVarDataByName(varName);
434 else
435 throw "Invalid variable name";
436 }
437
438 //
439 //
440 //
441 const QuadMaskInfo& FinleyElements::getQuadMask(int fsCode) const
442 {
443 if (fsCode == FINLEY_REDUCED_ELEMENTS ||
444 fsCode == FINLEY_REDUCED_FACE_ELEMENTS ||
445 fsCode == FINLEY_REDUCED_CONTACT_ELEMENTS_1) {
446 return reducedQuadMask;
447 } else {
448 return quadMask;
449 }
450 }
451
452 //
453 //
454 //
455 void FinleyElements::reorderArray(IntVec& v, const IntVec& idx,
456 int elementsPerIndex)
457 {
458 IntVec newArray(v.size());
459 IntVec::iterator arrIt = newArray.begin();
460 IntVec::const_iterator idxIt;
461 if (elementsPerIndex == 1) {
462 for (idxIt=idx.begin(); idxIt!=idx.end(); idxIt++) {
463 *arrIt++ = v[*idxIt];
464 }
465 } else {
466 for (idxIt=idx.begin(); idxIt!=idx.end(); idxIt++) {
467 int i = *idxIt;
468 copy(&v[i*elementsPerIndex], &v[(i+1)*elementsPerIndex], arrIt);
469 arrIt += elementsPerIndex;
470 }
471 }
472 v.swap(newArray);
473 }
474
475 //
476 //
477 //
478 void FinleyElements::buildReducedElements(const FinleyElementInfo& f)
479 {
480 // create the node list for the new element type
481 IntVec reducedNodes(f.reducedElementSize*numElements, 0);
482
483 IntVec::iterator reducedIt = reducedNodes.begin();
484 IntVec::const_iterator origIt;
485 for (origIt=nodes.begin(); origIt!=nodes.end();
486 origIt+=nodesPerElement)
487 {
488 copy(origIt, origIt+f.reducedElementSize, reducedIt);
489 reducedIt += f.reducedElementSize;
490 }
491
492 if (f.elementFactor > 1) {
493 // replace each element by multiple smaller ones which will be the
494 // new 'full' elements, whereas the original ones are the reduced
495 // elements, e.g.:
496 //
497 // new reduced: new full:
498 // _________ _________
499 // | | | | |
500 // | | > |____|____|
501 // | | > | | |
502 // |_________| |____|____|
503 //
504
505 // create the reduced elements which are basically a copy of the
506 // current elements
507 reducedElements = FinleyElements_ptr(new FinleyElements(
508 "Reduced"+name, originalMesh));
509 reducedElements->nodes = reducedNodes;
510 reducedElements->numElements = numElements;
511 reducedElements->type = f.reducedElementType;
512 reducedElements->nodesPerElement = f.reducedElementSize;
513 reducedElements->owner = owner;
514 reducedElements->color = color;
515 reducedElements->ID = ID;
516 reducedElements->tag = tag;
517
518 // now update full element data
519 IntVec fullNodes(f.elementSize*f.elementFactor*numElements);
520 IntVec::iterator cellIt = fullNodes.begin();
521
522 owner.clear();
523 color.clear();
524 ID.clear();
525 tag.clear();
526 for (int i=0; i < numElements; i++) {
527 owner.insert(owner.end(), f.elementFactor, reducedElements->owner[i]);
528 color.insert(color.end(), f.elementFactor, reducedElements->color[i]);
529 ID.insert(ID.end(), f.elementFactor, reducedElements->ID[i]);
530 tag.insert(tag.end(), f.elementFactor, reducedElements->tag[i]);
531 for (int j=0; j < f.elementFactor*f.elementSize; j++)
532 *cellIt++ = nodes[nodesPerElement*i+f.multiCellIndices[j]];
533 }
534
535 nodes.swap(fullNodes);
536 nodesPerElement = f.elementSize;
537 numElements *= f.elementFactor;
538
539 } else {
540 // we merely converted element types and don't need reduced elements
541 // so just replace node list and type
542 nodes.swap(reducedNodes);
543 nodesPerElement = f.reducedElementSize;
544 type = f.reducedElementType;
545 }
546 }
547
548 //
549 //
550 //
551 IntVec FinleyElements::prepareGhostIndices(int ownIndex)
552 {
553 IntVec indexArray;
554 numGhostElements = 0;
555
556 // move indices of "ghost zones" to the end to be able to reorder
557 // data accordingly
558 for (int i=0; i<numElements; i++) {
559 if (owner[i] == ownIndex)
560 indexArray.push_back(i);
561 }
562
563 for (int i=0; i<numElements; i++) {
564 if (owner[i] != ownIndex) {
565 numGhostElements++;
566 indexArray.push_back(i);
567 }
568 }
569 return indexArray;
570 }
571
572 //
573 //
574 //
575 void FinleyElements::reorderGhostZones(int ownIndex)
576 {
577 IntVec indexArray = prepareGhostIndices(ownIndex);
578
579 // move "ghost data" to the end of the arrays
580 if (numGhostElements > 0) {
581 reorderArray(nodes, indexArray, nodesPerElement);
582 reorderArray(owner, indexArray, 1);
583 reorderArray(color, indexArray, 1);
584 reorderArray(ID, indexArray, 1);
585 reorderArray(tag, indexArray, 1);
586 }
587
588 if (reducedElements)
589 reducedElements->reorderGhostZones(ownIndex);
590 }
591
592 //
593 //
594 //
595 void FinleyElements::removeGhostZones(int ownIndex)
596 {
597 reorderGhostZones(ownIndex);
598
599 if (numGhostElements > 0) {
600 numElements -= numGhostElements;
601 nodes.resize(numElements*nodesPerElement);
602 owner.resize(numElements);
603 color.resize(numElements);
604 ID.resize(numElements);
605 tag.resize(numElements);
606 numGhostElements = 0;
607 }
608
609 if (reducedElements)
610 reducedElements->removeGhostZones(ownIndex);
611 }
612
613 //
614 //
615 //
616 void FinleyElements::buildMeshes()
617 {
618 // build a new mesh containing only the required nodes
619 if (numElements > 0) {
620 if (nodeMesh && nodeMesh->getNumNodes() > 0) {
621 FinleyNodes_ptr newMesh(new FinleyNodes(nodeMesh, nodes, name));
622 nodeMesh.swap(newMesh);
623 } else {
624 nodeMesh.reset(new FinleyNodes(originalMesh, nodes, name));
625 }
626 #ifdef _DEBUG
627 cout << nodeMesh->getName() << " has " << nodeMesh->getNumNodes()
628 << " nodes and " << numElements << " elements" << endl;
629 #endif
630 }
631
632 if (reducedElements)
633 reducedElements->buildMeshes();
634 }
635
636 //
637 //
638 //
639 void FinleyElements::writeConnectivityVTK(ostream& os)
640 {
641 if (numElements > 0) {
642 const IntVec& gNI = nodeMesh->getGlobalNodeIndices();
643 IntVec::const_iterator it;
644 int count = 1;
645 for (it=nodes.begin(); it!=nodes.end(); it++, count++) {
646 os << gNI[*it];
647 if (count % nodesPerElement == 0)
648 os << endl;
649 else
650 os << " ";
651 }
652 }
653 }
654
655 #if USE_SILO
656 //
657 //
658 //
659 inline int toSiloElementType(int type)
660 {
661 switch (type) {
662 case ZONETYPE_BEAM: return DB_ZONETYPE_BEAM;
663 case ZONETYPE_HEX: return DB_ZONETYPE_HEX;
664 case ZONETYPE_POLYGON: return DB_ZONETYPE_POLYGON;
665 case ZONETYPE_QUAD: return DB_ZONETYPE_QUAD;
666 case ZONETYPE_TET: return DB_ZONETYPE_TET;
667 case ZONETYPE_TRIANGLE: return DB_ZONETYPE_TRIANGLE;
668 }
669 return 0;
670 }
671 #endif
672
673 //
674 //
675 //
676 bool FinleyElements::writeToSilo(DBfile* dbfile, const string& siloPath,
677 const StringVec& labels,
678 const StringVec& units)
679 {
680 #if USE_SILO
681 if (numElements == 0)
682 return true;
683
684 int ret;
685
686 if (siloPath != "") {
687 ret = DBSetDir(dbfile, siloPath.c_str());
688 if (ret != 0)
689 return false;
690 }
691
692 // write out the full mesh in any case
693 nodeMesh->setSiloPath(siloPath);
694 string siloMeshNameStr = nodeMesh->getFullSiloName();
695 const char* siloMeshName = siloMeshNameStr.c_str();
696 int arraylen = numElements * nodesPerElement;
697 int eltype = toSiloElementType(type);
698
699 string varName = name + string("_zones");
700 ret = DBPutZonelist2(dbfile, varName.c_str(), numElements,
701 nodeMesh->getNumDims(), &nodes[0], arraylen, 0, 0,
702 numGhostElements, &eltype, &nodesPerElement, &numElements, 1, NULL);
703 if (ret == 0) {
704 CoordArray& coordbase = const_cast<CoordArray&>(nodeMesh->getCoords());
705 DBoptlist* optList = NULL;
706 int nOpts = labels.size()+units.size();
707 if (nOpts>0) {
708 optList = DBMakeOptlist(nOpts);
709 if (labels.size()>0)
710 DBAddOption(optList, DBOPT_XLABEL, (void*)labels[0].c_str());
711 if (labels.size()>1)
712 DBAddOption(optList, DBOPT_YLABEL, (void*)labels[1].c_str());
713 if (labels.size()>2)
714 DBAddOption(optList, DBOPT_ZLABEL, (void*)labels[2].c_str());
715 if (units.size()>0)
716 DBAddOption(optList, DBOPT_XUNITS, (void*)units[0].c_str());
717 if (units.size()>1)
718 DBAddOption(optList, DBOPT_YUNITS, (void*)units[1].c_str());
719 if (units.size()>2)
720 DBAddOption(optList, DBOPT_ZUNITS, (void*)units[2].c_str());
721 }
722 ret = DBPutUcdmesh(dbfile, siloMeshName,
723 nodeMesh->getNumDims(), NULL, &coordbase[0],
724 nodeMesh->getNumNodes(), numElements, varName.c_str(),
725 /*"facelist"*/NULL, DB_FLOAT, optList);
726 if (optList) {
727 DBFreeOptlist(optList);
728 }
729 }
730
731 // Point mesh is useful for debugging
732 if (0) {
733 CoordArray& coordbase = const_cast<CoordArray&>(nodeMesh->getCoords());
734 DBPutPointmesh(dbfile, "/pointmesh", nodeMesh->getNumDims(),
735 &coordbase[0], nodeMesh->getNumNodes(), DB_FLOAT, NULL);
736 }
737
738 if (ret != 0)
739 return false;
740
741 // write out the element-centered variables
742 varName = name + string("_Color");
743 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
744 (float*)&color[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
745 NULL);
746 if (ret == 0) {
747 varName = name + string("_Id");
748 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
749 (float*)&ID[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
750 NULL);
751 }
752 if (ret == 0) {
753 varName = name + string("_Owner");
754 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
755 (float*)&owner[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
756 NULL);
757 }
758 if (ret == 0) {
759 varName = name + string("_Tag");
760 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
761 (float*)&tag[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
762 NULL);
763 }
764
765 if (reducedElements) {
766 reducedElements->writeToSilo(dbfile, siloPath, labels, units);
767 }
768
769 // "Elements" is a special case
770 if (name == "Elements") {
771 nodeMesh->writeToSilo(dbfile);
772 }
773
774 return (ret == 0);
775
776 #else // !USE_SILO
777 return false;
778 #endif
779 }
780
781 //
782 //
783 //
784 FinleyElementInfo FinleyElements::getDudleyTypeInfo(Dudley_ElementTypeId typeId)
785 {
786 FinleyElementInfo ret;
787 ret.multiCellIndices = NULL;
788 ret.elementFactor = 1;
789 ret.useQuadNodes = false;
790 ret.quadDim = 0;
791
792 switch (typeId) {
793 case Dudley_Line2Face://untested
794 case Dudley_Point1://untested
795 cerr << "WARNING: Dudley type " <<typeId<< " is untested!" << endl;
796 ret.elementSize = 1;
797 ret.elementType = ZONETYPE_POLYGON;
798 break;
799
800 case Dudley_Tri3Face://untested
801 cerr << "WARNING: Dudley type " <<typeId<< " is untested!" << endl;
802 case Dudley_Line2:
803 ret.elementSize = ret.reducedElementSize = 2;
804 ret.elementType = ret.reducedElementType = ZONETYPE_BEAM;
805 break;
806
807 case Dudley_Tet4Face://untested
808 cerr << "WARNING: Dudley type " <<typeId<< " is untested!" << endl;
809 case Dudley_Tri3:
810 ret.elementSize = ret.reducedElementSize = 3;
811 ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
812 break;
813
814 case Dudley_Tet4:
815 ret.elementSize = ret.reducedElementSize = 4;
816 ret.elementType = ret.reducedElementType = ZONETYPE_TET;
817 break;
818
819 default:
820 cerr << "WARNING: Unknown Dudley Type " << typeId << endl;
821 break;
822 }
823 return ret;
824 }
825
826 //
827 //
828 //
829 FinleyElementInfo FinleyElements::getFinleyTypeInfo(Finley_ElementTypeId typeId)
830 {
831 FinleyElementInfo ret;
832 ret.multiCellIndices = NULL;
833 ret.elementFactor = 1;
834 ret.useQuadNodes = false;
835 ret.quadDim = 0;
836
837 switch (typeId) {
838 case Finley_Point1_Contact://untested
839 case Finley_Line2Face_Contact://untested
840 case Finley_Line3Face_Contact://untested
841 case Finley_Line2Face://untested
842 case Finley_Line3Face://untested
843 case Finley_Point1://untested
844 cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
845 ret.elementSize = 1;
846 ret.elementType = ZONETYPE_POLYGON;
847 break;
848
849 case Finley_Tri3Face://untested
850 cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
851 case Finley_Tri3Face_Contact:
852 case Finley_Line2_Contact:
853 case Finley_Rec4Face_Contact:
854 case Finley_Rec4Face:
855 case Finley_Line2:
856 ret.elementSize = ret.reducedElementSize = 2;
857 ret.elementType = ret.reducedElementType = ZONETYPE_BEAM;
858 break;
859
860 case Finley_Line3Macro:
861 ret.useQuadNodes = true;
862 ret.quadDim = 1;
863 // fall through
864 case Finley_Line3:
865 ret.multiCellIndices = line3indices;
866 ret.elementFactor = 2;
867 // fall through
868 case Finley_Line3_Contact:
869 case Finley_Tri6Face_Contact://untested
870 case Finley_Rec8Face_Contact:
871 case Finley_Tri6Face://untested
872 case Finley_Rec8Face:
873 //VTK_QUADRATIC_EDGE
874 ret.elementSize = ret.reducedElementSize = 2;
875 ret.elementType = ret.reducedElementType = ZONETYPE_BEAM;
876 break;
877
878 case Finley_Tet4Face_Contact://untested
879 case Finley_Tet4Face://untested
880 cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
881 case Finley_Tri3_Contact:
882 case Finley_Tri3:
883 ret.elementSize = ret.reducedElementSize = 3;
884 ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
885 break;
886
887 case Finley_Rec4_Contact:
888 case Finley_Hex8Face_Contact:
889 case Finley_Hex8Face:
890 case Finley_Rec4:
891 ret.elementSize = ret.reducedElementSize = 4;
892 ret.elementType = ret.reducedElementType = ZONETYPE_QUAD;
893 break;
894
895 case Finley_Rec9:
896 case Finley_Rec9Macro:
897 ret.useQuadNodes = true;
898 ret.quadDim = 2;
899 ret.multiCellIndices = rec9indices;
900 ret.elementFactor = 4;
901 // fall through
902 case Finley_Rec9_Contact:
903 ret.elementSize = ret.reducedElementSize = 4;
904 ret.elementType = ret.reducedElementType = ZONETYPE_QUAD;
905 break;
906
907 case Finley_Tet4:
908 ret.elementSize = ret.reducedElementSize = 4;
909 ret.elementType = ret.reducedElementType = ZONETYPE_TET;
910 break;
911
912 case Finley_Tri6:
913 case Finley_Tri6Macro:
914 ret.useQuadNodes = true;
915 ret.quadDim = 2;
916 ret.multiCellIndices = tri6indices;
917 ret.elementFactor = 4;
918 // fall through
919 case Finley_Tri6_Contact:
920 case Finley_Tet10Face_Contact://untested
921 case Finley_Tet10Face://untested
922 //VTK_QUADRATIC_TRIANGLE
923 ret.elementSize = ret.reducedElementSize = 3;
924 ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
925 break;
926
927 case Finley_Rec8:
928 ret.multiCellIndices = rec8indices;
929 ret.elementFactor = 6;
930 // fall through
931 case Finley_Hex20Face:
932 case Finley_Rec8_Contact:
933 case Finley_Hex20Face_Contact:
934 //VTK_QUADRATIC_QUAD
935 ret.elementSize = 3;
936 ret.elementType = ZONETYPE_TRIANGLE;
937 ret.reducedElementSize = 4;
938 ret.reducedElementType = ZONETYPE_QUAD;
939 break;
940
941 case Finley_Tet10:
942 case Finley_Tet10Macro:
943 //VTK_QUADRATIC_TETRA
944 ret.useQuadNodes = true;
945 ret.quadDim = 3;
946 ret.multiCellIndices = tet10indices;
947 ret.elementFactor = 8;
948 ret.elementSize = ret.reducedElementSize = 4;
949 ret.elementType = ret.reducedElementType = ZONETYPE_TET;
950 break;
951
952 case Finley_Hex20:
953 //VTK_QUADRATIC_HEXAHEDRON
954 ret.multiCellIndices = hex20indices;
955 ret.elementFactor = 36;
956 ret.elementSize = 3;
957 ret.elementType = ZONETYPE_TRIANGLE;
958 ret.reducedElementSize = 8;
959 ret.reducedElementType = ZONETYPE_HEX;
960 break;
961
962 case Finley_Hex27:
963 case Finley_Hex27Macro:
964 ret.useQuadNodes = true;
965 ret.quadDim = 3;
966 ret.multiCellIndices = hex27indices;
967 ret.elementFactor = 8;
968 // fall through
969 case Finley_Hex8:
970 ret.elementSize = ret.reducedElementSize = 8;
971 ret.elementType = ret.reducedElementType = ZONETYPE_HEX;
972 break;
973
974 default:
975 cerr << "WARNING: Unknown Finley Type " << typeId << endl;
976 break;
977 }
978 return ret;
979 }
980
981 /////////////////////////////////
982 // Helpers for buildQuadMask() //
983 /////////////////////////////////
984
985 // returns true if |x-c| <= r, false otherwise
986 inline bool inside1D(float x, float c, float r)
987 {
988 return (ABS(x-c) <= r);
989 }
990
991 // returns true if |x-cx| <= r and |y-cy| <= r, false otherwise
992 inline bool inside2D(float x, float y, float cx, float cy, float r)
993 {
994 return (inside1D(x, cx, r) && inside1D(y, cy, r));
995 }
996
997 // returns true if |x-cx| <= r and |y-cy| <= r and |z-cz| <= r, false otherwise
998 inline bool inside3D(float x, float y, float z,
999 float cx, float cy, float cz, float r)
1000 {
1001 return (inside2D(x, y, cx, cy, r) && inside1D(z, cz, r));
1002 }
1003
1004 // returns true if d1 and d2 have the same sign or at least one of them is
1005 // close to 0, false otherwise
1006 inline bool sameSide(float d1, float d2)
1007 {
1008 const float TOL = 1.e-8f;
1009 return (ABS(d1) < TOL || ABS(d2) < TOL || d1*d2>=0.);
1010 }
1011
1012 // computes the determinant of the 4x4 matrix given by its elements m_ij
1013 static float det4x4(float m_00, float m_01, float m_02, float m_03,
1014 float m_10, float m_11, float m_12, float m_13,
1015 float m_20, float m_21, float m_22, float m_23,
1016 float m_30, float m_31, float m_32, float m_33)
1017 {
1018 float det1 = m_12 * m_23 - m_22 * m_13;
1019 float det2 = m_11 * m_23 - m_21 * m_13;
1020 float det3 = m_11 * m_22 - m_21 * m_12;
1021 float det4 = m_10 * m_23 - m_20 * m_13;
1022 float det5 = m_10 * m_22 - m_20 * m_12;
1023 float det6 = m_10 * m_21 - m_20 * m_11;
1024 return -m_30 * (m_01 * det1 - m_02 * det2 + m_03 * det3) +
1025 m_31 * (m_00 * det1 - m_02 * det4 + m_03 * det5) -
1026 m_32 * (m_00 * det2 - m_01 * det4 + m_03 * det6) +
1027 m_33 * (m_00 * det3 - m_01 * det5 + m_02 * det6);
1028 }
1029
1030 // returns true if point (x,y,z) is in or on the tetrahedron given by its
1031 // corner points p0, p1, p2 and p3, false otherwise.
1032 static bool pointInTet(float x, float y, float z,
1033 const float* p0, const float* p1,
1034 const float* p2, const float* p3)
1035 {
1036 float d0 = det4x4(
1037 p0[0], p0[1], p0[2], 1.f,
1038 p1[0], p1[1], p1[2], 1.f,
1039 p2[0], p2[1], p2[2], 1.f,
1040 p3[0], p3[1], p3[2], 1.f);
1041 float d1 = det4x4(
1042 x, y, z, 1.f,
1043 p1[0], p1[1], p1[2], 1.f,
1044 p2[0], p2[1], p2[2], 1.f,
1045 p3[0], p3[1], p3[2], 1.f);
1046 float d2 = det4x4(
1047 p0[0], p0[1], p0[2], 1.f,
1048 x, y, z, 1.f,
1049 p2[0], p2[1], p2[2], 1.f,
1050 p3[0], p3[1], p3[2], 1.f);
1051 float d3 = det4x4(
1052 p0[0], p0[1], p0[2], 1.f,
1053 p1[0], p1[1], p1[2], 1.f,
1054 x, y, z, 1.f,
1055 p3[0], p3[1], p3[2], 1.f);
1056 float d4 = det4x4(
1057 p0[0], p0[1], p0[2], 1.f,
1058 p1[0], p1[1], p1[2], 1.f,
1059 p2[0], p2[1], p2[2], 1.f,
1060 x, y, z, 1.f);
1061
1062 return (sameSide(d0,d1) && sameSide(d1,d2) &&
1063 sameSide(d2,d3) && sameSide(d3,d4));
1064 }
1065
1066 // returns true if point (x,y) is in or on the triangle given by its corner
1067 // points p0, p1 and p2, false otherwise
1068 static bool pointInTri(float x, float y,
1069 const float* p0, const float* p1, const float* p2)
1070 {
1071 const float TOL = 1.e-8f;
1072 float v0[2] = { p2[0]-p0[0], p2[1]-p0[1] };
1073 float v1[2] = { p1[0]-p0[0], p1[1]-p0[1] };
1074 float v2[2] = { x - p0[0], y - p0[1] };
1075
1076 float dot00 = v0[0]*v0[0]+v0[1]*v0[1];
1077 float dot01 = v0[0]*v1[0]+v0[1]*v1[1];
1078 float dot02 = v0[0]*v2[0]+v0[1]*v2[1];
1079 float dot11 = v1[0]*v1[0]+v1[1]*v1[1];
1080 float dot12 = v1[0]*v2[0]+v1[1]*v2[1];
1081 float invDenom = dot00*dot11 - dot01*dot01;
1082 if (ABS(invDenom) < TOL) invDenom = TOL;
1083 invDenom = 1./invDenom;
1084 float u = (dot11*dot02 - dot01*dot12) * invDenom;
1085 float v = (dot00*dot12 - dot01*dot02) * invDenom;
1086 return (u>=0.f) && (v>=0.f) && (u+v<=1.f);
1087 }
1088
1089 //
1090 //
1091 //
1092 QuadMaskInfo FinleyElements::buildQuadMask(const CoordArray& qnodes, int numQNodes)
1093 {
1094 QuadMaskInfo qmi;
1095 if (numQNodes == 0)
1096 return qmi;
1097
1098 if (finleyTypeId == Finley_Line3Macro) {
1099 for (int i=0; i<elementFactor; i++) {
1100 const float bounds[] = { .25, .75 };
1101 IntVec m(numQNodes, 0);
1102 int hits = 0;
1103 for (size_t j=0; j<numQNodes; j++) {
1104 if (inside1D(qnodes[0][j], bounds[i], .25)) {
1105 m[j] = 1;
1106 hits++;
1107 }
1108 }
1109 qmi.mask.push_back(m);
1110 if (hits == 0)
1111 qmi.factor.push_back(1);
1112 else
1113 qmi.factor.push_back(hits);
1114 }
1115 } else if ((finleyTypeId == Finley_Tri6) || (finleyTypeId == Finley_Tri6Macro)) {
1116 for (int i=0; i<elementFactor; i++) {
1117 const float bounds[][2] = { { 0., 0. }, { 1., 0. },
1118 { 0., 1. }, { .5, 0. },
1119 { .5, .5 }, { 0., .5 } };
1120 const size_t* nodeIdx = &tri6indices[i*nodesPerElement];
1121 IntVec m(numQNodes, 0);
1122 int hits = 0;
1123 for (size_t j=0; j<numQNodes; j++) {
1124 // check if point j is in triangle i
1125 if (pointInTri(qnodes[0][j], qnodes[1][j],
1126 bounds[nodeIdx[0]], bounds[nodeIdx[1]],
1127 bounds[nodeIdx[2]])) {
1128 m[j] = 1;
1129 hits++;
1130 }
1131 }
1132 if (hits == 0) {
1133 // if an element does not contain any quadrature points we
1134 // simply average over all data points within that element
1135 m = IntVec(numQNodes, 1);
1136 qmi.factor.push_back(numQNodes);
1137 } else {
1138 qmi.factor.push_back(hits);
1139 }
1140 qmi.mask.push_back(m);
1141 }
1142 } else if ((finleyTypeId == Finley_Tet10) || (finleyTypeId == Finley_Tet10Macro)) {
1143 for (int i=0; i<elementFactor; i++) {
1144 const float bounds[][3] = {
1145 { 0., 0., 0. },
1146 { 1., 0., 0. },
1147 { 0., 0., 1. },
1148 { 0., 1., 0. },
1149 { .5, 0., 0. },
1150 { .5, 0., .5 },
1151 { 0., 0., .5 },
1152 { 0., .5, 0. },
1153 { .5, .5, 0. },
1154 { 0., .5, .5 }
1155 };
1156 // need to reorder the elements
1157 const size_t elNumIdx[] = { 0,1,2,4,3,5,6,7 };
1158 const size_t* nodeIdx = &tet10indices[elNumIdx[i]*nodesPerElement];
1159 IntVec m(numQNodes, 0);
1160 int hits = 0;
1161 for (size_t j=0; j<numQNodes; j++) {
1162 // check if point j is in tetrahedron i
1163 if (pointInTet(qnodes[0][j], qnodes[1][j], qnodes[2][j],
1164 bounds[nodeIdx[0]], bounds[nodeIdx[1]],
1165 bounds[nodeIdx[2]], bounds[nodeIdx[3]])) {
1166 m[j] = 1;
1167 hits++;
1168 }
1169 }
1170 if (hits == 0) {
1171 // if an element does not contain any quadrature points we
1172 // simply average over all data points within that element
1173 m = IntVec(numQNodes, 1);
1174 qmi.factor.push_back(numQNodes);
1175 } else {
1176 qmi.factor.push_back(hits);
1177 }
1178 qmi.mask.push_back(m);
1179 }
1180 } else if ((finleyTypeId == Finley_Rec9) || (finleyTypeId == Finley_Rec9Macro)) {
1181 for (int i=0; i<elementFactor; i++) {
1182 const float bounds[][2] = { { 0.25, 0.25 }, { 0.75, 0.25 },
1183 { 0.25, 0.75 }, { 0.75, 0.75 } };
1184 IntVec m(numQNodes, 0);
1185 int hits = 0;
1186 for (size_t j=0; j<numQNodes; j++) {
1187 if (inside2D(qnodes[0][j], qnodes[1][j],
1188 bounds[i][0], bounds[i][1], .25)) {
1189 m[j] = 1;
1190 hits++;
1191 }
1192 }
1193 qmi.mask.push_back(m);
1194 if (hits == 0)
1195 qmi.factor.push_back(1);
1196 else
1197 qmi.factor.push_back(hits);
1198 }
1199 } else if ((finleyTypeId == Finley_Hex27) || (finleyTypeId == Finley_Hex27Macro) ){
1200 for (int i=0; i<elementFactor; i++) {
1201 const float bounds[][3] = {
1202 { 0.25, 0.25, 0.25 }, { 0.75, 0.25, 0.25 },
1203 { 0.25, 0.75, 0.25 }, { 0.75, 0.75, 0.25 },
1204 { 0.25, 0.25, 0.75 }, { 0.75, 0.25, 0.75 },
1205 { 0.25, 0.75, 0.75 }, { 0.75, 0.75, 0.75 } };
1206 IntVec m(numQNodes, 0);
1207 int hits = 0;
1208 for (size_t j=0; j<numQNodes; j++) {
1209 if (inside3D(qnodes[0][j], qnodes[1][j], qnodes[2][j],
1210 bounds[i][0], bounds[i][1], bounds[i][2], .25)) {
1211 m[j] = 1;
1212 hits++;
1213 }
1214 }
1215 qmi.mask.push_back(m);
1216 if (hits == 0)
1217 qmi.factor.push_back(1);
1218 else
1219 qmi.factor.push_back(hits);
1220 }
1221 }
1222 return qmi;
1223 }
1224
1225 } // namespace weipa
1226

  ViewVC Help
Powered by ViewVC 1.1.26