Commit d07ee80d authored by Kenneth Leiter's avatar Kenneth Leiter
Browse files

ENH: Add Hexahedron() to Hexahedron_125() Topology Conversion.

parent 4f65d6c6
......@@ -137,6 +137,12 @@ boost::shared_ptr<const XdmfTopologyType> XdmfTopologyType::Hexahedron_64()
return p;
}
boost::shared_ptr<const XdmfTopologyType> XdmfTopologyType::Hexahedron_125()
{
static boost::shared_ptr<const XdmfTopologyType> p(new XdmfTopologyType(125, "Hexahedron_125", Quartic));
return p;
}
boost::shared_ptr<const XdmfTopologyType> XdmfTopologyType::Mixed()
{
static boost::shared_ptr<const XdmfTopologyType> p(new XdmfTopologyType(0, "Mixed", Arbitrary));
......
......@@ -27,6 +27,10 @@
* Pyramid_13
* Wedge_15
* Hexahedron_20
* Hexahedron_24
* Hexahedron_27
* Hexahedron_64
* Hexahedron_125
* Mixed
* TwoDSMesh
* TwoDRectMesh
......@@ -44,7 +48,7 @@ public:
friend class XdmfTopology;
enum CellType {
NoCellType, Linear, Quadratic, Cubic, Arbitrary, Structured
NoCellType, Linear, Quadratic, Cubic, Quartic, Arbitrary, Structured
};
// Supported Xdmf Topology Types
......@@ -68,6 +72,7 @@ public:
static boost::shared_ptr<const XdmfTopologyType> Hexahedron_24();
static boost::shared_ptr<const XdmfTopologyType> Hexahedron_27();
static boost::shared_ptr<const XdmfTopologyType> Hexahedron_64();
static boost::shared_ptr<const XdmfTopologyType> Hexahedron_125();
static boost::shared_ptr<const XdmfTopologyType> Mixed();
static boost::shared_ptr<const XdmfTopologyType> TwoDSMesh();
static boost::shared_ptr<const XdmfTopologyType> TwoDRectMesh();
......
......@@ -10,82 +10,94 @@
#include "XdmfTopologyConverter.hpp"
#include "XdmfTopologyType.hpp"
XdmfTopologyConverter::XdmfTopologyConverter()
{
}
/**
* PIMPL
*/
class XdmfTopologyConverter::XdmfTopologyConverterImpl {
XdmfTopologyConverter::~XdmfTopologyConverter()
{
}
struct XdmfTopologyConverter::PointComparison {
bool operator()(const std::vector<double> & point1, const std::vector<double> & point2) const
{
double epsilon = 1e-6;
for(unsigned int i=0; i<3; ++i)
{
if(fabs(point1[i] - point2[i]) > epsilon)
{
return point1[i] < point2[i];
}
}
return false;
}
};
class XdmfTopologyConverter::HexahedronToHexahedron64 {
public:
HexahedronToHexahedron64()
XdmfTopologyConverterImpl()
{
}
inline void computeInteriorPoints(std::vector<double> & leftPoint, std::vector<double> & rightPoint, std::vector<double> & point1, std::vector<double> & point2) const
~XdmfTopologyConverterImpl()
{
leftPoint[0] = (1.0/3.0)*(point2[0] + 2*point1[0]);
leftPoint[1] = (1.0/3.0)*(point2[1] + 2*point1[1]);
leftPoint[2] = (1.0/3.0)*(point2[2] + 2*point1[2]);
}
rightPoint[0] = (1.0/3.0)*(2*point2[0] + point1[0]);
rightPoint[1] = (1.0/3.0)*(2*point2[1] + point1[1]);
rightPoint[2] = (1.0/3.0)*(2*point2[2] + point1[2]);
struct PointComparison {
bool operator()(const std::vector<double> & point1, const std::vector<double> & point2) const
{
double epsilon = 1e-6;
for(unsigned int i=0; i<3; ++i)
{
if(fabs(point1[i] - point2[i]) > epsilon)
{
return point1[i] < point2[i];
}
}
return false;
}
};
inline void insertPointWithoutCheck(const std::vector<double> & newPoint, const boost::shared_ptr<XdmfArray> & newConnectivity, const boost::shared_ptr<XdmfArray> & newPoints) const
{
newConnectivity->pushBack(newPoints->getSize() / 3);
newPoints->pushBack(newPoint[0]);
newPoints->pushBack(newPoint[1]);
newPoints->pushBack(newPoint[2]);
}
inline void computeLeftPoint(std::vector<double> & leftPoint, std::vector<double> & point1, std::vector<double> & point2) const
inline void insertPointWithCheck(const std::vector<double> & newPoint, std::map<std::vector<double>, unsigned int, PointComparison> & coordToIdMap, const boost::shared_ptr<XdmfArray> & newConnectivity, const boost::shared_ptr<XdmfArray> & newPoints) const
{
std::map<std::vector<double>, unsigned int>::const_iterator iter = coordToIdMap.find(newPoint);
if(iter == coordToIdMap.end())
{
// Not inserted before
coordToIdMap[newPoint] = newPoints->getSize() / 3;;
insertPointWithoutCheck(newPoint, newConnectivity, newPoints);
}
else
{
newConnectivity->pushBack(iter->second);
}
}
class HexahedronToHexahedron_64;
class HexahedronToHexahedron_125;
class Hexahedron_64ToHexahedron;
};
class XdmfTopologyConverter::XdmfTopologyConverterImpl::HexahedronToHexahedron_64 : public XdmfTopologyConverter::XdmfTopologyConverterImpl {
public:
HexahedronToHexahedron_64()
{
leftPoint[0] = (1.0/3.0)*(point2[0] + 2*point1[0]);
leftPoint[1] = (1.0/3.0)*(point2[1] + 2*point1[1]);
leftPoint[2] = (1.0/3.0)*(point2[2] + 2*point1[2]);
}
inline void computeRightPoint(std::vector<double> & rightPoint, std::vector<double> & point1, std::vector<double> & point2) const
inline void computeInteriorPoints(std::vector<double> & leftPoint, std::vector<double> & rightPoint, const std::vector<double> & point1, const std::vector<double> & point2) const
{
rightPoint[0] = (1.0/3.0)*(2*point2[0] + point1[0]);
rightPoint[1] = (1.0/3.0)*(2*point2[1] + point1[1]);
rightPoint[2] = (1.0/3.0)*(2*point2[2] + point1[2]);
leftPoint[0] = (1.0/3.0)*(point2[0] + 2*point1[0]);
leftPoint[1] = (1.0/3.0)*(point2[1] + 2*point1[1]);
leftPoint[2] = (1.0/3.0)*(point2[2] + 2*point1[2]);
rightPoint[0] = (1.0/3.0)*(2*point2[0] + point1[0]);
rightPoint[1] = (1.0/3.0)*(2*point2[1] + point1[1]);
rightPoint[2] = (1.0/3.0)*(2*point2[2] + point1[2]);
}
inline void insertPointWithoutCheck(std::vector<double> & newPoint, const boost::shared_ptr<XdmfArray> & newConnectivity, const boost::shared_ptr<XdmfArray> & newPoints) const
inline void computeLeftPoint(std::vector<double> & leftPoint, const std::vector<double> & point1, const std::vector<double> & point2) const
{
newConnectivity->pushBack(newPoints->getSize() / 3);
newPoints->pushBack(newPoint[0]);
newPoints->pushBack(newPoint[1]);
newPoints->pushBack(newPoint[2]);
leftPoint[0] = (1.0/3.0)*(point2[0] + 2*point1[0]);
leftPoint[1] = (1.0/3.0)*(point2[1] + 2*point1[1]);
leftPoint[2] = (1.0/3.0)*(point2[2] + 2*point1[2]);
}
inline void insertPointWithCheck(std::vector<double> & newPoint, std::map<std::vector<double>, unsigned int, PointComparison> & coordToIdMap, const boost::shared_ptr<XdmfArray> & newConnectivity, const boost::shared_ptr<XdmfArray> & newPoints) const
inline void computeRightPoint(std::vector<double> & rightPoint, const std::vector<double> & point1, const std::vector<double> & point2) const
{
std::map<std::vector<double>, unsigned int>::const_iterator iter = coordToIdMap.find(newPoint);
if(iter == coordToIdMap.end())
{
// Not inserted before
coordToIdMap[newPoint] = newPoints->getSize() / 3;;
insertPointWithoutCheck(newPoint, newConnectivity, newPoints);
}
else
{
newConnectivity->pushBack(iter->second);
}
rightPoint[0] = (1.0/3.0)*(2*point2[0] + point1[0]);
rightPoint[1] = (1.0/3.0)*(2*point2[1] + point1[1]);
rightPoint[2] = (1.0/3.0)*(2*point2[2] + point1[2]);
}
boost::shared_ptr<XdmfGrid> convert(const boost::shared_ptr<XdmfGrid> gridToConvert) const
......@@ -120,11 +132,11 @@ public:
std::vector<double> rightPoint(3);
std::map<std::vector<double>, unsigned int, PointComparison> coordToIdMap;
std::vector<std::vector<double> > localNodes(44, std::vector<double>(3));
for(unsigned int i=0; i<gridToConvert->getTopology()->getNumberElements(); ++i)
{
// Fill localNodes with original coordinate information.
std::vector<std::vector<double> > localNodes(8, std::vector<double>(3));
localNodes.reserve(44);
for(int j=0; j<8; ++j)
{
gridToConvert->getGeometry()->getArray()->getValuesCopy(gridToConvert->getTopology()->getArray()->getValueCopy<unsigned int>(8*i + j) * 3, &localNodes[j][0], 3);
......@@ -136,139 +148,139 @@ public:
// Case 0
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[0], localNodes[1]);
localNodes.push_back(leftPoint);
localNodes.push_back(rightPoint);
localNodes[8] = leftPoint;
localNodes[9] = rightPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 1
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[1], localNodes[2]);
localNodes.push_back(leftPoint);
localNodes.push_back(rightPoint);
localNodes[10] = leftPoint;
localNodes[11] = rightPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 2
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[2], localNodes[3]);
localNodes.push_back(leftPoint);
localNodes.push_back(rightPoint);
localNodes[12] = leftPoint;
localNodes[13] = rightPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 3
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[3], localNodes[0]);
localNodes.push_back(leftPoint);
localNodes.push_back(rightPoint);
localNodes[14] = leftPoint;
localNodes[15] = rightPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 4
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[4], localNodes[5]);
localNodes.push_back(leftPoint);
localNodes.push_back(rightPoint);
localNodes[16] = leftPoint;
localNodes[17] = rightPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 5
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[5], localNodes[6]);
localNodes.push_back(leftPoint);
localNodes.push_back(rightPoint);
localNodes[18] = leftPoint;
localNodes[19] = rightPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 6
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[6], localNodes[7]);
localNodes.push_back(leftPoint);
localNodes.push_back(rightPoint);
localNodes[20] = leftPoint;
localNodes[21] = rightPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 7
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[7], localNodes[4]);
localNodes.push_back(leftPoint);
localNodes.push_back(rightPoint);
localNodes[22] = leftPoint;
localNodes[23] = rightPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 8
computeLeftPoint(leftPoint, localNodes[0], localNodes[4]);
localNodes.push_back(leftPoint);
localNodes[24] = leftPoint;
insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
// Case 9
this->computeLeftPoint(leftPoint, localNodes[1], localNodes[5]);
localNodes.push_back(leftPoint);
localNodes[25] = leftPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
// Case 10
this->computeLeftPoint(leftPoint, localNodes[2], localNodes[6]);
localNodes.push_back(leftPoint);
localNodes[26] = leftPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
// Case 11
this->computeLeftPoint(leftPoint, localNodes[3], localNodes[7]);
localNodes.push_back(leftPoint);
localNodes[27] = leftPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
// Case 12
this->computeRightPoint(leftPoint, localNodes[0], localNodes[4]);
localNodes.push_back(leftPoint);
localNodes[28] = leftPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
// Case 13
this->computeRightPoint(leftPoint, localNodes[1], localNodes[5]);
localNodes.push_back(leftPoint);
localNodes[29] = leftPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
// Case 14
this->computeRightPoint(leftPoint, localNodes[2], localNodes[6]);
localNodes.push_back(leftPoint);
localNodes[30] = leftPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
// Case 15
this->computeRightPoint(leftPoint, localNodes[3], localNodes[7]);
localNodes.push_back(leftPoint);
localNodes[31] = leftPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
// Case 16
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[27], localNodes[24]);
localNodes.push_back(leftPoint);
localNodes.push_back(rightPoint);
localNodes[32] = leftPoint;
localNodes[33] = rightPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 17
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[25], localNodes[26]);
localNodes.push_back(leftPoint);
localNodes.push_back(rightPoint);
localNodes[34] = leftPoint;
localNodes[35] = rightPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 18
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[24], localNodes[25]);
localNodes.push_back(leftPoint);
localNodes.push_back(rightPoint);
localNodes[36] = leftPoint;
localNodes[37] = rightPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 19
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[26], localNodes[27]);
localNodes.push_back(leftPoint);
localNodes.push_back(rightPoint);
localNodes[38] = leftPoint;
localNodes[39] = rightPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 20
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[31], localNodes[28]);
localNodes.push_back(leftPoint);
localNodes.push_back(rightPoint);
localNodes[40] = leftPoint;
localNodes[41] = rightPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 21
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[29], localNodes[30]);
localNodes.push_back(leftPoint);
localNodes.push_back(rightPoint);
localNodes[42] = leftPoint;
localNodes[43] = rightPoint;
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
......@@ -287,49 +299,504 @@ public:
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 25
// Case 25
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[11], localNodes[14]);
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 26
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[23], localNodes[18]);
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 27
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[19], localNodes[22]);
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 28
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[33], localNodes[34]);
this->insertPointWithoutCheck(leftPoint, newConnectivity, newPoints);
this->insertPointWithoutCheck(rightPoint, newConnectivity, newPoints);
// Case 29
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[35], localNodes[32]);
this->insertPointWithoutCheck(leftPoint, newConnectivity, newPoints);
this->insertPointWithoutCheck(rightPoint, newConnectivity, newPoints);
// Case 30
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[41], localNodes[42]);
this->insertPointWithoutCheck(leftPoint, newConnectivity, newPoints);
this->insertPointWithoutCheck(rightPoint, newConnectivity, newPoints);
// Case 31
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[43], localNodes[40]);
this->insertPointWithoutCheck(leftPoint, newConnectivity, newPoints);
this->insertPointWithoutCheck(rightPoint, newConnectivity, newPoints);
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 26
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[23], localNodes[18]);
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 27
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[19], localNodes[22]);
this->insertPointWithCheck(leftPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(rightPoint, coordToIdMap, newConnectivity, newPoints);
// Case 28
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[33], localNodes[34]);
this->insertPointWithoutCheck(leftPoint, newConnectivity, newPoints);
this->insertPointWithoutCheck(rightPoint, newConnectivity, newPoints);
// Case 29
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[35], localNodes[32]);
this->insertPointWithoutCheck(leftPoint, newConnectivity, newPoints);
this->insertPointWithoutCheck(rightPoint, newConnectivity, newPoints);
// Case 30
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[41], localNodes[42]);
this->insertPointWithoutCheck(leftPoint, newConnectivity, newPoints);
this->insertPointWithoutCheck(rightPoint, newConnectivity, newPoints);
// Case 31
this->computeInteriorPoints(leftPoint, rightPoint, localNodes[43], localNodes[40]);
this->insertPointWithoutCheck(leftPoint, newConnectivity, newPoints);
this->insertPointWithoutCheck(rightPoint, newConnectivity, newPoints);
}
return toReturn;
}
};
class XdmfTopologyConverter::XdmfTopologyConverterImpl::HexahedronToHexahedron_125 : public XdmfTopologyConverter::XdmfTopologyConverterImpl {
public:
HexahedronToHexahedron_125()
{
}
inline void computeInteriorPoints(std::vector<double> & quarterPoint, std::vector<double> & midPoint, std::vector<double> & threeQuarterPoint, const std::vector<double> & point1, const std::vector<double> & point2) const
{
quarterPoint[0] = (1.0/4.0)*(point2[0] + 3*point1[0]);
quarterPoint[1] = (1.0/4.0)*(point2[1] + 3*point1[1]);
quarterPoint[2] = (1.0/4.0)*(point2[2] + 3*point1[2]);
midPoint[0] = (1.0/2.0)*(point2[0] + point1[0]);
midPoint[1] = (1.0/2.0)*(point2[1] + point1[1]);
midPoint[2] = (1.0/2.0)*(point2[2] + point1[2]);
threeQuarterPoint[0] = (1.0/4.0)*(3.0*point2[0] + point1[0]);
threeQuarterPoint[1] = (1.0/4.0)*(3.0*point2[1] + point1[1]);
threeQuarterPoint[2] = (1.0/4.0)*(3.0*point2[2] + point1[2]);
}
inline void computeQuarterPoint (std::vector<double> & quarterPoint, const std::vector<double> & point1, const std::vector<double> & point2) const
{
quarterPoint[0] = (1.0/4.0)*(point2[0] + 3*point1[0]);
quarterPoint[1] = (1.0/4.0)*(point2[1] + 3*point1[1]);
quarterPoint[2] = (1.0/4.0)*(point2[2] + 3*point1[2]);
}
inline void computeMidPoint(std::vector<double> & midPoint, const std::vector<double> & point1, const std::vector<double> & point2) const
{
midPoint[0] = (1.0/2.0)*(point2[0] + point1[0]);
midPoint[1] = (1.0/2.0)*(point2[1] + point1[1]);
midPoint[2] = (1.0/2.0)*(point2[2] + point1[2]);
}
inline void computeThreeQuarterPoint(std::vector<double> & threeQuarterPoint, const std::vector<double> & point1, const std::vector<double> & point2) const
{
threeQuarterPoint[0] = (1.0/4.0)*(3.0*point2[0] + point1[0]);
threeQuarterPoint[1] = (1.0/4.0)*(3.0*point2[1] + point1[1]);
threeQuarterPoint[2] = (1.0/4.0)*(3.0*point2[2] + point1[2]);
}
boost::shared_ptr<XdmfGrid> convert(const boost::shared_ptr<XdmfGrid> gridToConvert) const
{
boost::shared_ptr<XdmfGrid> toReturn = XdmfGrid::New();
toReturn->setName(gridToConvert->getName());
toReturn->getGeometry()->setType(gridToConvert->getGeometry()->getType());
toReturn->getTopology()->setType(XdmfTopologyType::Hexahedron_125());
boost::shared_ptr<XdmfArray> newPoints = toReturn->getGeometry()->getArray();
newPoints->initialize(gridToConvert->getGeometry()->getArray()->getType());
newPoints->resize(gridToConvert->getGeometry()->getArray()->getSize(), 0);
if(!gridToConvert->getGeometry()->getArray()->isInitialized())
{
gridToConvert->getGeometry()->getArray()->read();
}
// Copy all geometry values from old grid into new grid because we are keeping all old points.
newPoints->copyValues(0, gridToConvert->getGeometry()->getArray(), 0, gridToConvert->getGeometry()->getArray()->getSize());
boost::shared_ptr<XdmfArray> newConnectivity = toReturn->getTopology()->getArray();
newConnectivity->initialize(gridToConvert->getTopology()->getArray()->getType());
newConnectivity->reserve(125 * gridToConvert->getTopology()->getNumberElements());
if(!gridToConvert->getTopology()->getArray()->isInitialized())
{
gridToConvert->getTopology()->getArray()->read();
}
std::vector<double> quarterPoint(3);
std::vector<double> midPoint(3);
std::vector<double> threeQuarterPoint(3);
std::map<std::vector<double>, unsigned int, PointComparison> coordToIdMap;
std::vector<std::vector<double> > localNodes(80, std::vector<double>(3));
for(unsigned int i=0; i<gridToConvert->getTopology()->getNumberElements(); ++i)
{
// Fill localNodes with original coordinate information.
for(int j=0; j<8; ++j)
{
gridToConvert->getGeometry()->getArray()->getValuesCopy(gridToConvert->getTopology()->getArray()->getValueCopy<unsigned int>(8*i + j) * 3, &localNodes[j][0], 3);
}
// Add old connectivity information to newConnectivity.
newConnectivity->resize(newConnectivity->getSize() + 8, 0);
newConnectivity->copyValues(125*i, gridToConvert->getTopology()->getArray(), 8*i, 8);
// Case 0
this->computeInteriorPoints(quarterPoint, midPoint, threeQuarterPoint, localNodes[0], localNodes[1]);
localNodes[8] = quarterPoint;
localNodes[9] = midPoint;
localNodes[10] = threeQuarterPoint;
this->insertPointWithCheck(quarterPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(midPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(threeQuarterPoint, coordToIdMap, newConnectivity, newPoints);
// Case 1
this->computeInteriorPoints(quarterPoint, midPoint, threeQuarterPoint, localNodes[1], localNodes[2]);
localNodes[11] = quarterPoint;
localNodes[12] = midPoint;
localNodes[13] = threeQuarterPoint;
this->insertPointWithCheck(quarterPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(midPoint, coordToIdMap, newConnectivity, newPoints);
this->insertPointWithCheck(threeQuarterPoint, coordToIdMap, newConnectivity, newPoints);
// Case 2
this->computeInteriorPoints(quarterPoint, midPoint, threeQuarterPoint, localNodes[2], localNodes[3]);