/[escript]/branches/split/weipa/src/FinleyElements.cpp
ViewVC logotype

Contents of /branches/split/weipa/src/FinleyElements.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4724 - (show annotations)
Thu Mar 6 05:22:12 2014 UTC (5 years, 1 month ago) by jfenwick
File size: 39577 byte(s)
Work towards parallel domains

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

  ViewVC Help
Powered by ViewVC 1.1.26