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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4492 - (show annotations)
Tue Jul 2 01:44:11 2013 UTC (5 years, 9 months ago) by caltinay
File size: 39474 byte(s)
Finley changes that were held back while in release mode - moved more stuff
into finley namespace.

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

  ViewVC Help
Powered by ViewVC 1.1.26