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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (8 years, 8 months ago) by jfenwick
File size: 39054 byte(s)
Merging dudley and scons updates from branches

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_Contact://untested
850 case Finley_Tri3Face://untested
851 cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
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_Line3:
861 case Finley_Line3Macro:
862 ret.multiCellIndices = line3indices;
863 ret.elementFactor = 2;
864 // fall through
865 case Finley_Line3_Contact:
866 case Finley_Tri6Face_Contact://untested
867 case Finley_Rec8Face_Contact:
868 case Finley_Tri6Face://untested
869 case Finley_Rec8Face:
870 //VTK_QUADRATIC_EDGE
871 ret.elementSize = ret.reducedElementSize = 2;
872 ret.elementType = ret.reducedElementType = ZONETYPE_BEAM;
873 break;
874
875 case Finley_Tet4Face_Contact://untested
876 case Finley_Tet4Face://untested
877 cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
878 case Finley_Tri3_Contact:
879 case Finley_Tri3:
880 ret.elementSize = ret.reducedElementSize = 3;
881 ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
882 break;
883
884 case Finley_Rec4_Contact:
885 case Finley_Hex8Face_Contact:
886 case Finley_Hex8Face:
887 case Finley_Rec4:
888 ret.elementSize = ret.reducedElementSize = 4;
889 ret.elementType = ret.reducedElementType = ZONETYPE_QUAD;
890 break;
891
892 case Finley_Rec9:
893 case Finley_Rec9Macro:
894 ret.useQuadNodes = true;
895 ret.quadDim = 2;
896 ret.multiCellIndices = rec9indices;
897 ret.elementFactor = 4;
898 // fall through
899 case Finley_Rec9_Contact:
900 ret.elementSize = ret.reducedElementSize = 4;
901 ret.elementType = ret.reducedElementType = ZONETYPE_QUAD;
902 break;
903
904 case Finley_Tet4:
905 ret.elementSize = ret.reducedElementSize = 4;
906 ret.elementType = ret.reducedElementType = ZONETYPE_TET;
907 break;
908
909 case Finley_Tri6:
910 case Finley_Tri6Macro:
911 ret.useQuadNodes = true;
912 ret.quadDim = 2;
913 ret.multiCellIndices = tri6indices;
914 ret.elementFactor = 4;
915 // fall through
916 case Finley_Tri6_Contact:
917 case Finley_Tet10Face_Contact://untested
918 case Finley_Tet10Face://untested
919 //VTK_QUADRATIC_TRIANGLE
920 ret.elementSize = ret.reducedElementSize = 3;
921 ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
922 break;
923
924 case Finley_Rec8:
925 ret.multiCellIndices = rec8indices;
926 ret.elementFactor = 6;
927 // fall through
928 case Finley_Hex20Face:
929 case Finley_Rec8_Contact:
930 case Finley_Hex20Face_Contact:
931 //VTK_QUADRATIC_QUAD
932 ret.elementSize = 3;
933 ret.elementType = ZONETYPE_TRIANGLE;
934 ret.reducedElementSize = 4;
935 ret.reducedElementType = ZONETYPE_QUAD;
936 break;
937
938 case Finley_Tet10:
939 case Finley_Tet10Macro:
940 //VTK_QUADRATIC_TETRA
941 ret.useQuadNodes = true;
942 ret.quadDim = 3;
943 ret.multiCellIndices = tet10indices;
944 ret.elementFactor = 8;
945 ret.elementSize = ret.reducedElementSize = 4;
946 ret.elementType = ret.reducedElementType = ZONETYPE_TET;
947 break;
948
949 case Finley_Hex20:
950 //VTK_QUADRATIC_HEXAHEDRON
951 ret.multiCellIndices = hex20indices;
952 ret.elementFactor = 36;
953 ret.elementSize = 3;
954 ret.elementType = ZONETYPE_TRIANGLE;
955 ret.reducedElementSize = 8;
956 ret.reducedElementType = ZONETYPE_HEX;
957 break;
958
959 case Finley_Hex27:
960 case Finley_Hex27Macro:
961 ret.useQuadNodes = true;
962 ret.quadDim = 3;
963 ret.multiCellIndices = hex27indices;
964 ret.elementFactor = 8;
965 // fall through
966 case Finley_Hex8:
967 ret.elementSize = ret.reducedElementSize = 8;
968 ret.elementType = ret.reducedElementType = ZONETYPE_HEX;
969 break;
970
971 default:
972 cerr << "WARNING: Unknown Finley Type " << typeId << endl;
973 break;
974 }
975 return ret;
976 }
977
978 /////////////////////////////////
979 // Helpers for buildQuadMask() //
980 /////////////////////////////////
981
982 // returns true if |x-c| <= r, false otherwise
983 inline bool inside1D(float x, float c, float r)
984 {
985 return (ABS(x-c) <= r);
986 }
987
988 // returns true if |x-cx| <= r and |y-cy| <= r, false otherwise
989 inline bool inside2D(float x, float y, float cx, float cy, float r)
990 {
991 return (inside1D(x, cx, r) && inside1D(y, cy, r));
992 }
993
994 // returns true if |x-cx| <= r and |y-cy| <= r and |z-cz| <= r, false otherwise
995 inline bool inside3D(float x, float y, float z,
996 float cx, float cy, float cz, float r)
997 {
998 return (inside2D(x, y, cx, cy, r) && inside1D(z, cz, r));
999 }
1000
1001 // returns true if d1 and d2 have the same sign or at least one of them is
1002 // close to 0, false otherwise
1003 inline bool sameSide(float d1, float d2)
1004 {
1005 const float TOL = 1.e-8f;
1006 return (ABS(d1) < TOL || ABS(d2) < TOL || d1*d2>=0.);
1007 }
1008
1009 // computes the determinant of the 4x4 matrix given by its elements m_ij
1010 static float det4x4(float m_00, float m_01, float m_02, float m_03,
1011 float m_10, float m_11, float m_12, float m_13,
1012 float m_20, float m_21, float m_22, float m_23,
1013 float m_30, float m_31, float m_32, float m_33)
1014 {
1015 float det1 = m_12 * m_23 - m_22 * m_13;
1016 float det2 = m_11 * m_23 - m_21 * m_13;
1017 float det3 = m_11 * m_22 - m_21 * m_12;
1018 float det4 = m_10 * m_23 - m_20 * m_13;
1019 float det5 = m_10 * m_22 - m_20 * m_12;
1020 float det6 = m_10 * m_21 - m_20 * m_11;
1021 return -m_30 * (m_01 * det1 - m_02 * det2 + m_03 * det3) +
1022 m_31 * (m_00 * det1 - m_02 * det4 + m_03 * det5) -
1023 m_32 * (m_00 * det2 - m_01 * det4 + m_03 * det6) +
1024 m_33 * (m_00 * det3 - m_01 * det5 + m_02 * det6);
1025 }
1026
1027 // returns true if point (x,y,z) is in or on the tetrahedron given by its
1028 // corner points p0, p1, p2 and p3, false otherwise.
1029 static bool pointInTet(float x, float y, float z,
1030 const float* p0, const float* p1,
1031 const float* p2, const float* p3)
1032 {
1033 float d0 = det4x4(
1034 p0[0], p0[1], p0[2], 1.f,
1035 p1[0], p1[1], p1[2], 1.f,
1036 p2[0], p2[1], p2[2], 1.f,
1037 p3[0], p3[1], p3[2], 1.f);
1038 float d1 = det4x4(
1039 x, y, z, 1.f,
1040 p1[0], p1[1], p1[2], 1.f,
1041 p2[0], p2[1], p2[2], 1.f,
1042 p3[0], p3[1], p3[2], 1.f);
1043 float d2 = det4x4(
1044 p0[0], p0[1], p0[2], 1.f,
1045 x, y, z, 1.f,
1046 p2[0], p2[1], p2[2], 1.f,
1047 p3[0], p3[1], p3[2], 1.f);
1048 float d3 = det4x4(
1049 p0[0], p0[1], p0[2], 1.f,
1050 p1[0], p1[1], p1[2], 1.f,
1051 x, y, z, 1.f,
1052 p3[0], p3[1], p3[2], 1.f);
1053 float d4 = det4x4(
1054 p0[0], p0[1], p0[2], 1.f,
1055 p1[0], p1[1], p1[2], 1.f,
1056 p2[0], p2[1], p2[2], 1.f,
1057 x, y, z, 1.f);
1058
1059 return (sameSide(d0,d1) && sameSide(d1,d2) &&
1060 sameSide(d2,d3) && sameSide(d3,d4));
1061 }
1062
1063 // returns true if point (x,y) is in or on the triangle given by its corner
1064 // points p0, p1 and p2, false otherwise
1065 static bool pointInTri(float x, float y,
1066 const float* p0, const float* p1, const float* p2)
1067 {
1068 const float TOL = 1.e-8f;
1069 float v0[2] = { p2[0]-p0[0], p2[1]-p0[1] };
1070 float v1[2] = { p1[0]-p0[0], p1[1]-p0[1] };
1071 float v2[2] = { x - p0[0], y - p0[1] };
1072
1073 float dot00 = v0[0]*v0[0]+v0[1]*v0[1];
1074 float dot01 = v0[0]*v1[0]+v0[1]*v1[1];
1075 float dot02 = v0[0]*v2[0]+v0[1]*v2[1];
1076 float dot11 = v1[0]*v1[0]+v1[1]*v1[1];
1077 float dot12 = v1[0]*v2[0]+v1[1]*v2[1];
1078 float invDenom = dot00*dot11 - dot01*dot01;
1079 if (ABS(invDenom) < TOL) invDenom = TOL;
1080 invDenom = 1./invDenom;
1081 float u = (dot11*dot02 - dot01*dot12) * invDenom;
1082 float v = (dot00*dot12 - dot01*dot02) * invDenom;
1083 return (u>=0.f) && (v>=0.f) && (u+v<=1.f);
1084 }
1085
1086 //
1087 //
1088 //
1089 QuadMaskInfo FinleyElements::buildQuadMask(const CoordArray& qnodes, int numQNodes)
1090 {
1091 QuadMaskInfo qmi;
1092 if (numQNodes == 0)
1093 return qmi;
1094
1095 if (finleyTypeId == Finley_Line3Macro) {
1096 for (int i=0; i<elementFactor; i++) {
1097 const float bounds[] = { 0.25, 0.75 };
1098 IntVec m(numQNodes, 0);
1099 int hits = 0;
1100 for (size_t j=0; j<numQNodes; j++) {
1101 if (inside1D(qnodes[0][j], bounds[i], .25)) {
1102 m[j] = 1;
1103 hits++;
1104 }
1105 }
1106 qmi.mask.push_back(m);
1107 if (hits == 0)
1108 qmi.factor.push_back(1);
1109 else
1110 qmi.factor.push_back(hits);
1111 }
1112 } else if ((finleyTypeId == Finley_Tri6) || (finleyTypeId == Finley_Tri6Macro)) {
1113 for (int i=0; i<elementFactor; i++) {
1114 const float bounds[][2] = { { 0., 0. }, { 1., 0. },
1115 { 0., 1. }, { .5, 0. },
1116 { .5, .5 }, { 0., .5 } };
1117 const size_t* nodeIdx = &tri6indices[i*nodesPerElement];
1118 IntVec m(numQNodes, 0);
1119 int hits = 0;
1120 for (size_t j=0; j<numQNodes; j++) {
1121 // check if point j is in triangle i
1122 if (pointInTri(qnodes[0][j], qnodes[1][j],
1123 bounds[nodeIdx[0]], bounds[nodeIdx[1]],
1124 bounds[nodeIdx[2]])) {
1125 m[j] = 1;
1126 hits++;
1127 }
1128 }
1129 if (hits == 0) {
1130 // if an element does not contain any quadrature points we
1131 // simply average over all data points within that element
1132 m = IntVec(numQNodes, 1);
1133 qmi.factor.push_back(numQNodes);
1134 qmi.factor.push_back(1);
1135 } else {
1136 qmi.factor.push_back(hits);
1137 }
1138 qmi.mask.push_back(m);
1139 }
1140 } else if ((finleyTypeId == Finley_Tet10) || (finleyTypeId == Finley_Tet10Macro)) {
1141 for (int i=0; i<elementFactor; i++) {
1142 const float bounds[][3] = {
1143 { 0., 0., 0. },
1144 { 1., 0., 0. },
1145 { 0., 0., 1. },
1146 { 0., 1., 0. },
1147 { .5, 0., 0. },
1148 { .5, 0., .5 },
1149 { 0., 0., .5 },
1150 { 0., .5, 0. },
1151 { .5, .5, 0. },
1152 { 0., .5, .5 }
1153 };
1154 // need to reorder the elements
1155 const size_t elNumIdx[] = { 0,1,2,4,3,5,6,7 };
1156 const size_t* nodeIdx = &tet10indices[elNumIdx[i]*nodesPerElement];
1157 IntVec m(numQNodes, 0);
1158 int hits = 0;
1159 for (size_t j=0; j<numQNodes; j++) {
1160 // check if point j is in tetrahedron i
1161 if (pointInTet(qnodes[0][j], qnodes[1][j], qnodes[2][j],
1162 bounds[nodeIdx[0]], bounds[nodeIdx[1]],
1163 bounds[nodeIdx[2]], bounds[nodeIdx[3]])) {
1164 m[j] = 1;
1165 hits++;
1166 }
1167 }
1168 if (hits == 0) {
1169 // if an element does not contain any quadrature points we
1170 // simply average over all data points within that element
1171 m = IntVec(numQNodes, 1);
1172 qmi.factor.push_back(numQNodes);
1173 } else {
1174 qmi.factor.push_back(hits);
1175 }
1176 qmi.mask.push_back(m);
1177 }
1178 } else if ((finleyTypeId == Finley_Rec9) || (finleyTypeId == Finley_Rec9Macro)) {
1179 for (int i=0; i<elementFactor; i++) {
1180 const float bounds[][2] = { { 0.25, 0.25 }, { 0.75, 0.25 },
1181 { 0.25, 0.75 }, { 0.75, 0.75 } };
1182 IntVec m(numQNodes, 0);
1183 int hits = 0;
1184 for (size_t j=0; j<numQNodes; j++) {
1185 if (inside2D(qnodes[0][j], qnodes[1][j],
1186 bounds[i][0], bounds[i][1], 0.25)) {
1187 m[j] = 1;
1188 hits++;
1189 }
1190 }
1191 qmi.mask.push_back(m);
1192 if (hits == 0)
1193 qmi.factor.push_back(1);
1194 else
1195 qmi.factor.push_back(hits);
1196 }
1197 } else if ((finleyTypeId == Finley_Hex27) || (finleyTypeId == Finley_Hex27Macro) ){
1198 for (int i=0; i<elementFactor; i++) {
1199 const float bounds[][3] = {
1200 { 0.25, 0.25, 0.25 }, { 0.75, 0.25, 0.25 },
1201 { 0.25, 0.75, 0.25 }, { 0.75, 0.75, 0.25 },
1202 { 0.25, 0.25, 0.75 }, { 0.75, 0.25, 0.75 },
1203 { 0.25, 0.75, 0.75 }, { 0.75, 0.75, 0.75 } };
1204 IntVec m(numQNodes, 0);
1205 int hits = 0;
1206 for (size_t j=0; j<numQNodes; j++) {
1207 if (inside3D(qnodes[0][j], qnodes[1][j], qnodes[2][j],
1208 bounds[i][0], bounds[i][1], bounds[i][2], 0.25)) {
1209 m[j] = 1;
1210 hits++;
1211 }
1212 }
1213 qmi.mask.push_back(m);
1214 if (hits == 0)
1215 qmi.factor.push_back(1);
1216 else
1217 qmi.factor.push_back(hits);
1218 }
1219 }
1220 return qmi;
1221 }
1222
1223 } // namespace weipa
1224

  ViewVC Help
Powered by ViewVC 1.1.26