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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3792 - (show annotations)
Wed Feb 1 06:16:25 2012 UTC (6 years, 11 months ago) by caltinay
File size: 39293 byte(s)
Merged ripley rectangular domain into trunk.

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 #elif not defined(ABS)
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, bool writeMeshData)
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 if enabled
742 if (writeMeshData) {
743 varName = name + string("_Color");
744 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
745 (float*)&color[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
746 NULL);
747 if (ret == 0) {
748 varName = name + string("_Id");
749 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
750 (float*)&ID[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
751 NULL);
752 }
753 if (ret == 0) {
754 varName = name + string("_Owner");
755 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
756 (float*)&owner[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
757 NULL);
758 }
759 if (ret == 0) {
760 varName = name + string("_Tag");
761 ret = DBPutUcdvar1(dbfile, varName.c_str(), siloMeshName,
762 (float*)&tag[0], numElements, NULL, 0, DB_INT, DB_ZONECENT,
763 NULL);
764 }
765 }
766
767 if (reducedElements) {
768 reducedElements->writeToSilo(dbfile, siloPath, labels, units, writeMeshData);
769 }
770
771 // "Elements" is a special case
772 if (writeMeshData && name == "Elements") {
773 nodeMesh->writeToSilo(dbfile);
774 }
775
776 return (ret == 0);
777
778 #else // !USE_SILO
779 return false;
780 #endif
781 }
782
783 //
784 //
785 //
786 FinleyElementInfo FinleyElements::getDudleyTypeInfo(Dudley_ElementTypeId typeId)
787 {
788 FinleyElementInfo ret;
789 ret.multiCellIndices = NULL;
790 ret.elementFactor = 1;
791 ret.useQuadNodes = false;
792 ret.quadDim = 0;
793
794 switch (typeId) {
795 case Dudley_Line2Face://untested
796 case Dudley_Point1://untested
797 cerr << "WARNING: Dudley type " <<typeId<< " is untested!" << endl;
798 ret.elementSize = 1;
799 ret.elementType = ZONETYPE_POLYGON;
800 break;
801
802 case Dudley_Tri3Face://untested
803 cerr << "WARNING: Dudley type " <<typeId<< " is untested!" << endl;
804 case Dudley_Line2:
805 ret.elementSize = ret.reducedElementSize = 2;
806 ret.elementType = ret.reducedElementType = ZONETYPE_BEAM;
807 break;
808
809 case Dudley_Tet4Face://untested
810 cerr << "WARNING: Dudley type " <<typeId<< " is untested!" << endl;
811 case Dudley_Tri3:
812 ret.elementSize = ret.reducedElementSize = 3;
813 ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
814 break;
815
816 case Dudley_Tet4:
817 ret.elementSize = ret.reducedElementSize = 4;
818 ret.elementType = ret.reducedElementType = ZONETYPE_TET;
819 break;
820
821 default:
822 cerr << "WARNING: Unknown Dudley Type " << typeId << endl;
823 break;
824 }
825 return ret;
826 }
827
828 //
829 //
830 //
831 FinleyElementInfo FinleyElements::getFinleyTypeInfo(Finley_ElementTypeId typeId)
832 {
833 FinleyElementInfo ret;
834 ret.multiCellIndices = NULL;
835 ret.elementFactor = 1;
836 ret.useQuadNodes = false;
837 ret.quadDim = 0;
838
839 switch (typeId) {
840 case Finley_Point1_Contact://untested
841 case Finley_Line2Face_Contact://untested
842 case Finley_Line3Face_Contact://untested
843 case Finley_Line2Face://untested
844 case Finley_Line3Face://untested
845 case Finley_Point1://untested
846 cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
847 ret.elementSize = 1;
848 ret.elementType = ZONETYPE_POLYGON;
849 break;
850
851 case Finley_Tri3Face://untested
852 cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
853 case Finley_Tri3Face_Contact:
854 case Finley_Line2_Contact:
855 case Finley_Rec4Face_Contact:
856 case Finley_Rec4Face:
857 case Finley_Line2:
858 ret.elementSize = ret.reducedElementSize = 2;
859 ret.elementType = ret.reducedElementType = ZONETYPE_BEAM;
860 break;
861
862 case Finley_Line3Macro:
863 ret.useQuadNodes = true;
864 ret.quadDim = 1;
865 // fall through
866 case Finley_Line3:
867 ret.multiCellIndices = line3indices;
868 ret.elementFactor = 2;
869 // fall through
870 case Finley_Line3_Contact:
871 case Finley_Tri6Face_Contact://untested
872 case Finley_Rec8Face_Contact:
873 case Finley_Tri6Face://untested
874 case Finley_Rec8Face:
875 //VTK_QUADRATIC_EDGE
876 ret.elementSize = ret.reducedElementSize = 2;
877 ret.elementType = ret.reducedElementType = ZONETYPE_BEAM;
878 break;
879
880 case Finley_Tet4Face_Contact://untested
881 case Finley_Tet4Face://untested
882 cerr << "WARNING: Finley type " <<typeId<< " is untested!" << endl;
883 case Finley_Tri3_Contact:
884 case Finley_Tri3:
885 ret.elementSize = ret.reducedElementSize = 3;
886 ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
887 break;
888
889 case Finley_Rec4_Contact:
890 case Finley_Hex8Face_Contact:
891 case Finley_Hex8Face:
892 case Finley_Rec4:
893 ret.elementSize = ret.reducedElementSize = 4;
894 ret.elementType = ret.reducedElementType = ZONETYPE_QUAD;
895 break;
896
897 case Finley_Rec9:
898 case Finley_Rec9Macro:
899 ret.useQuadNodes = true;
900 ret.quadDim = 2;
901 ret.multiCellIndices = rec9indices;
902 ret.elementFactor = 4;
903 // fall through
904 case Finley_Rec9_Contact:
905 ret.elementSize = ret.reducedElementSize = 4;
906 ret.elementType = ret.reducedElementType = ZONETYPE_QUAD;
907 break;
908
909 case Finley_Tet4:
910 ret.elementSize = ret.reducedElementSize = 4;
911 ret.elementType = ret.reducedElementType = ZONETYPE_TET;
912 break;
913
914 case Finley_Tri6:
915 case Finley_Tri6Macro:
916 ret.useQuadNodes = true;
917 ret.quadDim = 2;
918 ret.multiCellIndices = tri6indices;
919 ret.elementFactor = 4;
920 // fall through
921 case Finley_Tri6_Contact:
922 case Finley_Tet10Face_Contact://untested
923 case Finley_Tet10Face://untested
924 //VTK_QUADRATIC_TRIANGLE
925 ret.elementSize = ret.reducedElementSize = 3;
926 ret.elementType = ret.reducedElementType = ZONETYPE_TRIANGLE;
927 break;
928
929 case Finley_Rec8:
930 ret.multiCellIndices = rec8indices;
931 ret.elementFactor = 6;
932 // fall through
933 case Finley_Hex20Face:
934 case Finley_Rec8_Contact:
935 case Finley_Hex20Face_Contact:
936 //VTK_QUADRATIC_QUAD
937 ret.elementSize = 3;
938 ret.elementType = ZONETYPE_TRIANGLE;
939 ret.reducedElementSize = 4;
940 ret.reducedElementType = ZONETYPE_QUAD;
941 break;
942
943 case Finley_Tet10:
944 case Finley_Tet10Macro:
945 //VTK_QUADRATIC_TETRA
946 ret.useQuadNodes = true;
947 ret.quadDim = 3;
948 ret.multiCellIndices = tet10indices;
949 ret.elementFactor = 8;
950 ret.elementSize = ret.reducedElementSize = 4;
951 ret.elementType = ret.reducedElementType = ZONETYPE_TET;
952 break;
953
954 case Finley_Hex20:
955 //VTK_QUADRATIC_HEXAHEDRON
956 ret.multiCellIndices = hex20indices;
957 ret.elementFactor = 36;
958 ret.elementSize = 3;
959 ret.elementType = ZONETYPE_TRIANGLE;
960 ret.reducedElementSize = 8;
961 ret.reducedElementType = ZONETYPE_HEX;
962 break;
963
964 case Finley_Hex27:
965 case Finley_Hex27Macro:
966 ret.useQuadNodes = true;
967 ret.quadDim = 3;
968 ret.multiCellIndices = hex27indices;
969 ret.elementFactor = 8;
970 // fall through
971 case Finley_Hex8:
972 ret.elementSize = ret.reducedElementSize = 8;
973 ret.elementType = ret.reducedElementType = ZONETYPE_HEX;
974 break;
975
976 default:
977 cerr << "WARNING: Unknown Finley Type " << typeId << endl;
978 break;
979 }
980 return ret;
981 }
982
983 /////////////////////////////////
984 // Helpers for buildQuadMask() //
985 /////////////////////////////////
986
987 // returns true if |x-c| <= r, false otherwise
988 inline bool inside1D(float x, float c, float r)
989 {
990 return (ABS(x-c) <= r);
991 }
992
993 // returns true if |x-cx| <= r and |y-cy| <= r, false otherwise
994 inline bool inside2D(float x, float y, float cx, float cy, float r)
995 {
996 return (inside1D(x, cx, r) && inside1D(y, cy, r));
997 }
998
999 // returns true if |x-cx| <= r and |y-cy| <= r and |z-cz| <= r, false otherwise
1000 inline bool inside3D(float x, float y, float z,
1001 float cx, float cy, float cz, float r)
1002 {
1003 return (inside2D(x, y, cx, cy, r) && inside1D(z, cz, r));
1004 }
1005
1006 // returns true if d1 and d2 have the same sign or at least one of them is
1007 // close to 0, false otherwise
1008 inline bool sameSide(float d1, float d2)
1009 {
1010 const float TOL = 1.e-8f;
1011 return (ABS(d1) < TOL || ABS(d2) < TOL || d1*d2>=0.);
1012 }
1013
1014 // computes the determinant of the 4x4 matrix given by its elements m_ij
1015 static float det4x4(float m_00, float m_01, float m_02, float m_03,
1016 float m_10, float m_11, float m_12, float m_13,
1017 float m_20, float m_21, float m_22, float m_23,
1018 float m_30, float m_31, float m_32, float m_33)
1019 {
1020 float det1 = m_12 * m_23 - m_22 * m_13;
1021 float det2 = m_11 * m_23 - m_21 * m_13;
1022 float det3 = m_11 * m_22 - m_21 * m_12;
1023 float det4 = m_10 * m_23 - m_20 * m_13;
1024 float det5 = m_10 * m_22 - m_20 * m_12;
1025 float det6 = m_10 * m_21 - m_20 * m_11;
1026 return -m_30 * (m_01 * det1 - m_02 * det2 + m_03 * det3) +
1027 m_31 * (m_00 * det1 - m_02 * det4 + m_03 * det5) -
1028 m_32 * (m_00 * det2 - m_01 * det4 + m_03 * det6) +
1029 m_33 * (m_00 * det3 - m_01 * det5 + m_02 * det6);
1030 }
1031
1032 // returns true if point (x,y,z) is in or on the tetrahedron given by its
1033 // corner points p0, p1, p2 and p3, false otherwise.
1034 static bool pointInTet(float x, float y, float z,
1035 const float* p0, const float* p1,
1036 const float* p2, const float* p3)
1037 {
1038 float d0 = det4x4(
1039 p0[0], p0[1], p0[2], 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 d1 = det4x4(
1044 x, y, z, 1.f,
1045 p1[0], p1[1], p1[2], 1.f,
1046 p2[0], p2[1], p2[2], 1.f,
1047 p3[0], p3[1], p3[2], 1.f);
1048 float d2 = det4x4(
1049 p0[0], p0[1], p0[2], 1.f,
1050 x, y, z, 1.f,
1051 p2[0], p2[1], p2[2], 1.f,
1052 p3[0], p3[1], p3[2], 1.f);
1053 float d3 = det4x4(
1054 p0[0], p0[1], p0[2], 1.f,
1055 p1[0], p1[1], p1[2], 1.f,
1056 x, y, z, 1.f,
1057 p3[0], p3[1], p3[2], 1.f);
1058 float d4 = det4x4(
1059 p0[0], p0[1], p0[2], 1.f,
1060 p1[0], p1[1], p1[2], 1.f,
1061 p2[0], p2[1], p2[2], 1.f,
1062 x, y, z, 1.f);
1063
1064 return (sameSide(d0,d1) && sameSide(d1,d2) &&
1065 sameSide(d2,d3) && sameSide(d3,d4));
1066 }
1067
1068 // returns true if point (x,y) is in or on the triangle given by its corner
1069 // points p0, p1 and p2, false otherwise
1070 static bool pointInTri(float x, float y,
1071 const float* p0, const float* p1, const float* p2)
1072 {
1073 const float TOL = 1.e-8f;
1074 float v0[2] = { p2[0]-p0[0], p2[1]-p0[1] };
1075 float v1[2] = { p1[0]-p0[0], p1[1]-p0[1] };
1076 float v2[2] = { x - p0[0], y - p0[1] };
1077
1078 float dot00 = v0[0]*v0[0]+v0[1]*v0[1];
1079 float dot01 = v0[0]*v1[0]+v0[1]*v1[1];
1080 float dot02 = v0[0]*v2[0]+v0[1]*v2[1];
1081 float dot11 = v1[0]*v1[0]+v1[1]*v1[1];
1082 float dot12 = v1[0]*v2[0]+v1[1]*v2[1];
1083 float invDenom = dot00*dot11 - dot01*dot01;
1084 if (ABS(invDenom) < TOL) invDenom = TOL;
1085 invDenom = 1.f/invDenom;
1086 float u = (dot11*dot02 - dot01*dot12) * invDenom;
1087 float v = (dot00*dot12 - dot01*dot02) * invDenom;
1088 return (u>=0.f) && (v>=0.f) && (u+v<=1.f);
1089 }
1090
1091 //
1092 //
1093 //
1094 QuadMaskInfo FinleyElements::buildQuadMask(const CoordArray& qnodes, int numQNodes)
1095 {
1096 QuadMaskInfo qmi;
1097 if (numQNodes == 0)
1098 return qmi;
1099
1100 if (finleyTypeId == Finley_Line3Macro) {
1101 for (int i=0; i<elementFactor; i++) {
1102 const float bounds[] = { .25, .75 };
1103 IntVec m(numQNodes, 0);
1104 int hits = 0;
1105 for (size_t j=0; j<numQNodes; j++) {
1106 if (inside1D(qnodes[0][j], bounds[i], .25)) {
1107 m[j] = 1;
1108 hits++;
1109 }
1110 }
1111 qmi.mask.push_back(m);
1112 if (hits == 0)
1113 qmi.factor.push_back(1);
1114 else
1115 qmi.factor.push_back(hits);
1116 }
1117 } else if ((finleyTypeId == Finley_Tri6) || (finleyTypeId == Finley_Tri6Macro)) {
1118 for (int i=0; i<elementFactor; i++) {
1119 const float bounds[][2] = { { 0., 0. }, { 1., 0. },
1120 { 0., 1. }, { .5, 0. },
1121 { .5, .5 }, { 0., .5 } };
1122 const size_t* nodeIdx = &tri6indices[i*nodesPerElement];
1123 IntVec m(numQNodes, 0);
1124 int hits = 0;
1125 for (size_t j=0; j<numQNodes; j++) {
1126 // check if point j is in triangle i
1127 if (pointInTri(qnodes[0][j], qnodes[1][j],
1128 bounds[nodeIdx[0]], bounds[nodeIdx[1]],
1129 bounds[nodeIdx[2]])) {
1130 m[j] = 1;
1131 hits++;
1132 }
1133 }
1134 if (hits == 0) {
1135 // if an element does not contain any quadrature points we
1136 // simply average over all data points within that element
1137 m = IntVec(numQNodes, 1);
1138 qmi.factor.push_back(numQNodes);
1139 } else {
1140 qmi.factor.push_back(hits);
1141 }
1142 qmi.mask.push_back(m);
1143 }
1144 } else if ((finleyTypeId == Finley_Tet10) || (finleyTypeId == Finley_Tet10Macro)) {
1145 for (int i=0; i<elementFactor; i++) {
1146 const float bounds[][3] = {
1147 { 0., 0., 0. },
1148 { 1., 0., 0. },
1149 { 0., 0., 1. },
1150 { 0., 1., 0. },
1151 { .5, 0., 0. },
1152 { .5, 0., .5 },
1153 { 0., 0., .5 },
1154 { 0., .5, 0. },
1155 { .5, .5, 0. },
1156 { 0., .5, .5 }
1157 };
1158 // need to reorder the elements
1159 const size_t elNumIdx[] = { 0,1,2,4,3,5,6,7 };
1160 const size_t* nodeIdx = &tet10indices[elNumIdx[i]*nodesPerElement];
1161 IntVec m(numQNodes, 0);
1162 int hits = 0;
1163 for (size_t j=0; j<numQNodes; j++) {
1164 // check if point j is in tetrahedron i
1165 if (pointInTet(qnodes[0][j], qnodes[1][j], qnodes[2][j],
1166 bounds[nodeIdx[0]], bounds[nodeIdx[1]],
1167 bounds[nodeIdx[2]], bounds[nodeIdx[3]])) {
1168 m[j] = 1;
1169 hits++;
1170 }
1171 }
1172 if (hits == 0) {
1173 // if an element does not contain any quadrature points we
1174 // simply average over all data points within that element
1175 m = IntVec(numQNodes, 1);
1176 qmi.factor.push_back(numQNodes);
1177 } else {
1178 qmi.factor.push_back(hits);
1179 }
1180 qmi.mask.push_back(m);
1181 }
1182 } else if ((finleyTypeId == Finley_Rec9) || (finleyTypeId == Finley_Rec9Macro)) {
1183 for (int i=0; i<elementFactor; i++) {
1184 const float bounds[][2] = { { 0.25, 0.25 }, { 0.75, 0.25 },
1185 { 0.25, 0.75 }, { 0.75, 0.75 } };
1186 IntVec m(numQNodes, 0);
1187 int hits = 0;
1188 for (size_t j=0; j<numQNodes; j++) {
1189 if (inside2D(qnodes[0][j], qnodes[1][j],
1190 bounds[i][0], bounds[i][1], .25)) {
1191 m[j] = 1;
1192 hits++;
1193 }
1194 }
1195 qmi.mask.push_back(m);
1196 if (hits == 0)
1197 qmi.factor.push_back(1);
1198 else
1199 qmi.factor.push_back(hits);
1200 }
1201 } else if ((finleyTypeId == Finley_Hex27) || (finleyTypeId == Finley_Hex27Macro) ){
1202 for (int i=0; i<elementFactor; i++) {
1203 const float bounds[][3] = {
1204 { 0.25, 0.25, 0.25 }, { 0.75, 0.25, 0.25 },
1205 { 0.25, 0.75, 0.25 }, { 0.75, 0.75, 0.25 },
1206 { 0.25, 0.25, 0.75 }, { 0.75, 0.25, 0.75 },
1207 { 0.25, 0.75, 0.75 }, { 0.75, 0.75, 0.75 } };
1208 IntVec m(numQNodes, 0);
1209 int hits = 0;
1210 for (size_t j=0; j<numQNodes; j++) {
1211 if (inside3D(qnodes[0][j], qnodes[1][j], qnodes[2][j],
1212 bounds[i][0], bounds[i][1], bounds[i][2], .25)) {
1213 m[j] = 1;
1214 hits++;
1215 }
1216 }
1217 qmi.mask.push_back(m);
1218 if (hits == 0)
1219 qmi.factor.push_back(1);
1220 else
1221 qmi.factor.push_back(hits);
1222 }
1223 }
1224 return qmi;
1225 }
1226
1227 } // namespace weipa
1228

  ViewVC Help
Powered by ViewVC 1.1.26