From 45377509670f2400bfe4ff93bbd2d128bef038e8 Mon Sep 17 00:00:00 2001 From: Kenneth Leiter Date: Thu, 8 Mar 2012 09:49:47 -0500 Subject: [PATCH] ENH: Add ability to convert Hexahedron to Hexahedron_27 in topology converter. --- utils/XdmfTopologyConverter.cpp | 429 ++++++++++++++++++++++++++++---- utils/XdmfTopologyConverter.hpp | 1 + 2 files changed, 380 insertions(+), 50 deletions(-) diff --git a/utils/XdmfTopologyConverter.cpp b/utils/XdmfTopologyConverter.cpp index 0a67f342..1fded2a5 100644 --- a/utils/XdmfTopologyConverter.cpp +++ b/utils/XdmfTopologyConverter.cpp @@ -41,6 +41,55 @@ // namespace { + void handleSetConversion(const shared_ptr gridToConvert, + const shared_ptr toReturn, + const std::map & oldIdToNewId, + const shared_ptr heavyDataWriter) + { + + for(unsigned int i=0; igetNumberSets(); ++i) { + const shared_ptr set = gridToConvert->getSet(i); + const shared_ptr setType = set->getType(); + if(setType == XdmfSetType::Cell()) { + toReturn->insert(set); + } + else if(setType == XdmfSetType::Node()) { + bool releaseSet = false; + if(!set->isInitialized()) { + set->read(); + releaseSet = true; + } + shared_ptr toReturnSet = XdmfSet::New(); + toReturnSet->setName(set->getName()); + toReturnSet->setType(set->getType()); + toReturnSet->initialize(set->getArrayType(), + set->getSize()); + + for(unsigned int i=0; igetSize(); ++i) { + const unsigned int nodeId = set->getValue(i); + std::map::const_iterator iter = + oldIdToNewId.find(nodeId); + if(iter == oldIdToNewId.end()) { + XdmfError::message(XdmfError::FATAL, + "Error converting hex node id set to hex_27 " + "node id set."); + } + toReturnSet->insert(i, iter->second); + } + if(releaseSet) { + set->release(); + } + + toReturn->insert(toReturnSet); + + if(heavyDataWriter) { + toReturnSet->accept(heavyDataWriter); + toReturnSet->release(); + } + } + } + } + // Classes that perform topology conversions. Converter is the root // base class. Tessellator is a subclass of Converter that deals // with cases where the mesh only needs to be tessellated to carry @@ -224,6 +273,327 @@ namespace { }; + class HexahedronToHexahedron27 : public Converter { + + public: + + HexahedronToHexahedron27() + { + } + + virtual ~HexahedronToHexahedron27() + { + } + + void + calculateMidPoint(std::vector & result, + const std::vector & point1, + const std::vector & point2) const + { + for (int i=0; i<3; i++) + result[i] = (point1[i]+point2[i]) / 2.0; + } + + shared_ptr + convert(const shared_ptr gridToConvert, + const shared_ptr topologyType, + const shared_ptr heavyDataWriter) const + { + + shared_ptr toReturn = XdmfUnstructuredGrid::New(); + toReturn->setName(gridToConvert->getName()); + + shared_ptr geometry = gridToConvert->getGeometry(); + shared_ptr toReturnGeometry = toReturn->getGeometry(); + + toReturnGeometry->setType(geometry->getType()); + toReturnGeometry->initialize(geometry->getArrayType()); + + bool releaseGeometry = false; + if(!geometry->isInitialized()) { + geometry->read(); + releaseGeometry = true; + } + + shared_ptr topology = gridToConvert->getTopology(); + shared_ptr toReturnTopology = toReturn->getTopology(); + + toReturn->getTopology()->setType(topologyType); + toReturnTopology->initialize(topology->getArrayType()); + toReturnTopology->reserve(topologyType->getNodesPerElement() * + topology->getNumberElements()); + + bool releaseTopology = false; + if(!topology->isInitialized()) { + topology->read(); + releaseTopology = true; + } + + std::map, unsigned int, PointComparison> coordToIdMap; + std::map oldIdToNewId; + + // allocate storage for values used in loop + unsigned int zeroIndex; + unsigned int oneIndex; + unsigned int twoIndex; + unsigned int threeIndex; + unsigned int fourIndex; + unsigned int fiveIndex; + unsigned int sixIndex; + unsigned int sevenIndex; + std::vector elementCorner0(3); + std::vector elementCorner1(3); + std::vector elementCorner2(3); + std::vector elementCorner3(3); + std::vector elementCorner4(3); + std::vector elementCorner5(3); + std::vector elementCorner6(3); + std::vector elementCorner7(3); + std::vector newPoint(3); + + unsigned int offset = 0; + unsigned int newIndex = 0; + for(unsigned int elem = 0; elemgetNumberElements(); ++elem) { + + // get indices of coner vertices of the element + zeroIndex = topology->getValue(offset++); + oneIndex = topology->getValue(offset++); + twoIndex = topology->getValue(offset++); + threeIndex = topology->getValue(offset++); + fourIndex = topology->getValue(offset++); + fiveIndex = topology->getValue(offset++); + sixIndex = topology->getValue(offset++); + sevenIndex = topology->getValue(offset++); + + // get locations of corner vertices of the element + geometry->getValues(zeroIndex * 3, + &(elementCorner0[0]), + 3); + geometry->getValues(oneIndex * 3, + &(elementCorner1[0]), + 3); + geometry->getValues(twoIndex * 3, + &(elementCorner2[0]), + 3); + geometry->getValues(threeIndex * 3, + &(elementCorner3[0]), + 3); + geometry->getValues(fourIndex * 3, + &(elementCorner4[0]), + 3); + geometry->getValues(fiveIndex * 3, + &(elementCorner5[0]), + 3); + geometry->getValues(sixIndex * 3, + &(elementCorner6[0]), + 3); + geometry->getValues(sevenIndex * 3, + &(elementCorner7[0]), + 3); + + // insert corner points + newIndex = this->insertPointWithCheck(elementCorner0, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + oldIdToNewId[zeroIndex] = newIndex; + newIndex = this->insertPointWithCheck(elementCorner1, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + oldIdToNewId[oneIndex] = newIndex; + newIndex = this->insertPointWithCheck(elementCorner2, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + oldIdToNewId[twoIndex] = newIndex; + newIndex = this->insertPointWithCheck(elementCorner3, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + oldIdToNewId[threeIndex] = newIndex; + newIndex = this->insertPointWithCheck(elementCorner4, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + oldIdToNewId[fourIndex] = newIndex; + newIndex = this->insertPointWithCheck(elementCorner5, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + oldIdToNewId[fiveIndex] = newIndex; + newIndex = this->insertPointWithCheck(elementCorner6, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + oldIdToNewId[sixIndex] = newIndex; + newIndex = this->insertPointWithCheck(elementCorner7, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + oldIdToNewId[sevenIndex] = newIndex; + + // insert additional points + calculateMidPoint(newPoint, + elementCorner0, + elementCorner1); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner1, + elementCorner2); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner2, + elementCorner3); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner0, + elementCorner3); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner4, + elementCorner5); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner5, + elementCorner6); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner6, + elementCorner7); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner4, + elementCorner7); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner0, + elementCorner4); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner1, + elementCorner5); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner2, + elementCorner6); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner3, + elementCorner7); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner0, + elementCorner7); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner1, + elementCorner6); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner1, + elementCorner4); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner2, + elementCorner7); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner0, + elementCorner2); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner4, + elementCorner6); + this->insertPointWithCheck(newPoint, + coordToIdMap, + toReturnTopology, + toReturnGeometry); + calculateMidPoint(newPoint, + elementCorner0, + elementCorner6); + this->insertPointWithoutCheck(newPoint, + toReturnTopology, + toReturnGeometry); + + } + + if(releaseTopology) { + topology->release(); + } + if(releaseGeometry) { + geometry->release(); + } + + if(heavyDataWriter) { + toReturnTopology->accept(heavyDataWriter); + toReturnTopology->release(); + toReturnGeometry->accept(heavyDataWriter); + toReturnGeometry->release(); + } + + handleSetConversion(gridToConvert, + toReturn, + oldIdToNewId, + heavyDataWriter); + + return toReturn; + + } + + }; + template class HexahedronToHighOrderHexahedron : public Converter { @@ -249,7 +619,6 @@ namespace { result[i] = point1[i]+scalar*(point2[i]-point1[i]); } - shared_ptr convert(const shared_ptr gridToConvert, const shared_ptr topologyType, @@ -276,7 +645,7 @@ namespace { toReturn->getTopology()->setType(topologyType); toReturnTopology->initialize(topology->getArrayType()); - toReturnTopology->reserve(mNodesPerElement * + toReturnTopology->reserve(topologyType->getNodesPerElement() * topology->getNumberElements()); bool releaseTopology = false; @@ -316,9 +685,7 @@ namespace { unsigned int offset = 0; for(unsigned int elem = 0; elemgetNumberElements(); ++elem) { - // // get indices of coner vertices of the element - // zeroIndex = topology->getValue(offset++); oneIndex = topology->getValue(offset++); twoIndex = topology->getValue(offset++); @@ -402,7 +769,7 @@ namespace { if((i == 0 || i == ORDER) || (j == 0 || j == ORDER) || (k == 0 || k == ORDER)) { - unsigned int newIndex = + const unsigned int newIndex = this->insertPointWithCheck(point, coordToIdMap, toReturnTopology, @@ -468,55 +835,16 @@ namespace { toReturnGeometry->release(); } - // handle sets - for(unsigned int i=0; igetNumberSets(); ++i) { - const shared_ptr set = gridToConvert->getSet(i); - const shared_ptr setType = set->getType(); - if(setType == XdmfSetType::Cell()) { - toReturn->insert(set); - } - else if(setType == XdmfSetType::Node()) { - bool releaseSet = false; - if(!set->isInitialized()) { - set->read(); - releaseSet = true; - } - shared_ptr toReturnSet = XdmfSet::New(); - toReturnSet->setName(set->getName()); - toReturnSet->setType(set->getType()); - toReturnSet->initialize(set->getArrayType(), - set->getSize()); - - for(unsigned int i=0; igetSize(); ++i) { - const unsigned int nodeId = set->getValue(i); - std::map::const_iterator iter = - oldIdToNewId.find(nodeId); - if(iter == oldIdToNewId.end()) { - XdmfError::message(XdmfError::FATAL, - "Error converting hex node id set to high " - "order node id set."); - } - toReturnSet->insert(i, iter->second); - } - if(releaseSet) { - set->release(); - } - - toReturn->insert(toReturnSet); + handleSetConversion(gridToConvert, + toReturn, + oldIdToNewId, + heavyDataWriter); - if(heavyDataWriter) { - toReturnSet->accept(heavyDataWriter); - toReturnSet->release(); - } - } - } return toReturn; } private: static const unsigned int mNodesPerEdge = ORDER + 1; - static const unsigned int mNodesPerElement = mNodesPerEdge * - mNodesPerEdge * mNodesPerEdge; static const double points[]; }; @@ -760,8 +1088,6 @@ namespace { private: static const unsigned int mNodesPerEdge = (ORDER + 1); - static const unsigned int mNodesPerElement = mNodesPerEdge * - mNodesPerEdge * mNodesPerEdge; }; @@ -809,6 +1135,9 @@ XdmfTopologyConverter::convert(const shared_ptr gridToConv Converter * converter = NULL; if(topologyTypeToConvert == XdmfTopologyType::Hexahedron()) { + if(topologyType == XdmfTopologyType::Hexahedron_27()) { + converter = new HexahedronToHexahedron27(); + } if(topologyType == XdmfTopologyType::Hexahedron_64()) { if(options == 1) { converter = new HexahedronToHighOrderHexahedron<3, true>(); diff --git a/utils/XdmfTopologyConverter.hpp b/utils/XdmfTopologyConverter.hpp index ef615bce..27b27039 100644 --- a/utils/XdmfTopologyConverter.hpp +++ b/utils/XdmfTopologyConverter.hpp @@ -47,6 +47,7 @@ class XdmfUnstructuredGrid; * form the new topology, no additional points are added. * * Currently supported conversions: + * Hexahedron to Hexahedron_27 * Hexahedron to Hexahedron_64 * Hexahedron to Hexahedron_125 * Hexahedron to Hexahedron_216 -- GitLab