Commit 5d773b0e authored by Kenneth Leiter's avatar Kenneth Leiter
Browse files

ENH: Add additional topology types.

Modify XdmfTopologyConverter and XdmfPartitioner to work with new
topology types.
parent 3ada9fd4
......@@ -220,29 +220,66 @@ XdmfTopologyType::Hexahedron_64()
}
shared_ptr<const XdmfTopologyType>
XdmfTopologyType::Hexahedron_64_GLL()
XdmfTopologyType::Hexahedron_125()
{
static shared_ptr<const XdmfTopologyType>
p(new XdmfTopologyType(64, "Hexahedron_64_GLL", Cubic, 0x35));
p(new XdmfTopologyType(125, "Hexahedron_125", Quartic, 0x34));
return p;
}
shared_ptr<const XdmfTopologyType>
XdmfTopologyType::Hexahedron_125()
XdmfTopologyType::Hexahedron_216()
{
static shared_ptr<const XdmfTopologyType>
p(new XdmfTopologyType(125, "Hexahedron_125", Quartic, 0x34));
p(new XdmfTopologyType(216, "Hexahedron_216", Quintic, 0x35));
return p;
}
shared_ptr<const XdmfTopologyType>
XdmfTopologyType::Hexahedron_125_GLL()
XdmfTopologyType::Hexahedron_343()
{
static shared_ptr<const XdmfTopologyType>
p(new XdmfTopologyType(125, "Hexahedron_125_GLL", Quartic, 0x36));
p(new XdmfTopologyType(343, "Hexahedron_343", Sextic, 0x36));
return p;
}
shared_ptr<const XdmfTopologyType>
XdmfTopologyType::Hexahedron_512()
{
static shared_ptr<const XdmfTopologyType>
p(new XdmfTopologyType(512, "Hexahedron_512", Septic, 0x37));
return p;
}
shared_ptr<const XdmfTopologyType>
XdmfTopologyType::Hexahedron_729()
{
static shared_ptr<const XdmfTopologyType>
p(new XdmfTopologyType(729, "Hexahedron_729", Octic, 0x38));
return p;
}
shared_ptr<const XdmfTopologyType>
XdmfTopologyType::Hexahedron_1000()
{
static shared_ptr<const XdmfTopologyType>
p(new XdmfTopologyType(1000, "Hexahedron_1000", Nonic, 0x39));
return p;
}
shared_ptr<const XdmfTopologyType>
XdmfTopologyType::Hexahedron_1331()
{
static shared_ptr<const XdmfTopologyType>
p(new XdmfTopologyType(1331, "Hexahedron_1331", Decic, 0x40));
return p;
}
shared_ptr<const XdmfTopologyType>
XdmfTopologyType::Mixed()
{
......@@ -320,14 +357,26 @@ XdmfTopologyType::New(const unsigned int id)
else if(id == XdmfTopologyType::Hexahedron_64()->getID()) {
return XdmfTopologyType::Hexahedron_64();
}
else if(id == XdmfTopologyType::Hexahedron_64_GLL()->getID()) {
return XdmfTopologyType::Hexahedron_64_GLL();
}
else if(id == XdmfTopologyType::Hexahedron_125()->getID()) {
return XdmfTopologyType::Hexahedron_125();
}
else if(id == XdmfTopologyType::Hexahedron_125_GLL()->getID()) {
return XdmfTopologyType::Hexahedron_125_GLL();
else if(id == XdmfTopologyType::Hexahedron_216()->getID()) {
return XdmfTopologyType::Hexahedron_216();
}
else if(id == XdmfTopologyType::Hexahedron_343()->getID()) {
return XdmfTopologyType::Hexahedron_343();
}
else if(id == XdmfTopologyType::Hexahedron_512()->getID()) {
return XdmfTopologyType::Hexahedron_512();
}
else if(id == XdmfTopologyType::Hexahedron_729()->getID()) {
return XdmfTopologyType::Hexahedron_729();
}
else if(id == XdmfTopologyType::Hexahedron_1000()->getID()) {
return XdmfTopologyType::Hexahedron_1000();
}
else if(id == XdmfTopologyType::Hexahedron_1331()->getID()) {
return XdmfTopologyType::Hexahedron_1331();
}
else if(id == XdmfTopologyType::Mixed()->getID()) {
return XdmfTopologyType::Mixed();
......@@ -440,6 +489,24 @@ XdmfTopologyType::New(const std::map<std::string, std::string> & itemProperties)
else if(typeVal.compare("HEXAHEDRON_125") == 0) {
return Hexahedron_125();
}
else if(typeVal.compare("HEXAHEDRON_216") == 0) {
return Hexahedron_216();
}
else if(typeVal.compare("HEXAHEDRON_343") == 0) {
return Hexahedron_343();
}
else if(typeVal.compare("HEXAHEDRON_512") == 0) {
return Hexahedron_512();
}
else if(typeVal.compare("HEXAHEDRON_729") == 0) {
return Hexahedron_729();
}
else if(typeVal.compare("HEXAHEDRON_1000") == 0) {
return Hexahedron_1000();
}
else if(typeVal.compare("HEXAHEDRON_1331") == 0) {
return Hexahedron_1331();
}
else if(typeVal.compare("MIXED") == 0) {
return Mixed();
}
......
......@@ -59,9 +59,13 @@
* Hexahedron_24 - 24 Node Bi-Quadratic Hexahedron
* Hexahedron_27 - 27 Node Tri-Quadratic Hexahedron
* Hexahedron_64 - 64 Node Tri-Cubic Hexahedron
* Hexahedron_64_GLL - 64 Node Spectral Tri-Cubic Hexahedron with Gauss-Lobatto-Legendre points.
* Hexahedron_125 - 125 Node Tri-Quartic Hexahedron
* Hexahedron_125_GLL - 125 Node Spectral Tri-Quartic Hexahedron with Gauss-Lobatto-Legendre points.
* Hexahedron_216 - 216 Node Tri-Quintic Hexahedron
* Hexahedron_343 - 343 Node Tri-Hexic Hexahedron
* Hexahedron_512 - 512 Node Tri-Septic Hexahedron
* Hexahedron_729 - 729 Node Tri-Octic Hexahedron
* Hexahedron_1000 - 1000 Node Tri-Nonic Hexahedron
* Hexahedron_1331 - 1331 Node Tri-Decic Hexahedron
* Mixed - Mixture of Unstructured Topologies
*/
class XDMF_EXPORT XdmfTopologyType : public XdmfItemProperty {
......@@ -73,13 +77,19 @@ public:
friend class XdmfTopology;
enum CellType {
NoCellType,
Linear,
Quadratic,
Cubic,
Quartic,
Arbitrary,
Structured
NoCellType = 0,
Linear = 1,
Quadratic = 2,
Cubic = 3,
Quartic = 4,
Quintic = 5,
Sextic = 6,
Septic = 7,
Octic = 8,
Nonic = 9,
Decic = 10,
Arbitrary = 100,
Structured = 101
};
/**
......@@ -109,9 +119,13 @@ public:
static shared_ptr<const XdmfTopologyType> Hexahedron_24();
static shared_ptr<const XdmfTopologyType> Hexahedron_27();
static shared_ptr<const XdmfTopologyType> Hexahedron_64();
static shared_ptr<const XdmfTopologyType> Hexahedron_64_GLL();
static shared_ptr<const XdmfTopologyType> Hexahedron_125();
static shared_ptr<const XdmfTopologyType> Hexahedron_125_GLL();
static shared_ptr<const XdmfTopologyType> Hexahedron_216();
static shared_ptr<const XdmfTopologyType> Hexahedron_343();
static shared_ptr<const XdmfTopologyType> Hexahedron_512();
static shared_ptr<const XdmfTopologyType> Hexahedron_729();
static shared_ptr<const XdmfTopologyType> Hexahedron_1000();
static shared_ptr<const XdmfTopologyType> Hexahedron_1331();
static shared_ptr<const XdmfTopologyType> Mixed();
/**
......
......@@ -113,7 +113,13 @@ XdmfPartitioner::partition(const shared_ptr<XdmfUnstructuredGrid> gridToPartitio
topologyType == XdmfTopologyType::Hexahedron_24() ||
topologyType == XdmfTopologyType::Hexahedron_27() ||
topologyType == XdmfTopologyType::Hexahedron_64() ||
topologyType == XdmfTopologyType::Hexahedron_125()) {
topologyType == XdmfTopologyType::Hexahedron_125() ||
topologyType == XdmfTopologyType::Hexahedron_216() ||
topologyType == XdmfTopologyType::Hexahedron_343() ||
topologyType == XdmfTopologyType::Hexahedron_512() ||
topologyType == XdmfTopologyType::Hexahedron_729() ||
topologyType == XdmfTopologyType::Hexahedron_1000() ||
topologyType == XdmfTopologyType::Hexahedron_1331()) {
metisElementType = 3;
nodesPerElement = 8;
}
......
......@@ -28,6 +28,8 @@
#include "XdmfGeometry.hpp"
#include "XdmfGeometryType.hpp"
#include "XdmfHeavyDataWriter.hpp"
#include "XdmfSet.hpp"
#include "XdmfSetType.hpp"
#include "XdmfTopology.hpp"
#include "XdmfTopologyConverter.hpp"
#include "XdmfTopologyType.hpp"
......@@ -39,19 +41,10 @@
//
namespace {
// 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 out the conversion
// (e.g. Hexahedron_64 to Hexahedron.
class Converter;
class Tessellator;
class HexahedronToHexahedron_64;
class HexahedronToHexahedron_64_GLL;
class HexahedronToHexahedron_125;
class HexahedronToHexahedron_125_GLL;
class Hexahedron_64ToHexahedron;
class Hexahedron_125ToHexahedron;
// 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
// out the conversion (e.g. Hexahedron_64 to Hexahedron).
class Converter {
......@@ -67,16 +60,19 @@ namespace {
virtual shared_ptr<XdmfUnstructuredGrid>
convert(const shared_ptr<XdmfUnstructuredGrid> gridToConvert,
const shared_ptr<const XdmfTopologyType> topologyType,
const shared_ptr<XdmfHeavyDataWriter> heavyDataWriter) const = 0;
protected:
struct PointComparison {
static const double epsilon = 1e-6;
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];
......@@ -86,18 +82,20 @@ namespace {
}
};
void
unsigned int
insertPointWithoutCheck(const std::vector<double> & newPoint,
const shared_ptr<XdmfTopology> & newConnectivity,
const shared_ptr<XdmfGeometry> & newPoints) const
{
newConnectivity->pushBack<unsigned int>(newPoints->getSize() / 3);
const unsigned int index = newPoints->getSize() / 3;
newConnectivity->pushBack(index);
newPoints->pushBack(newPoint[0]);
newPoints->pushBack(newPoint[1]);
newPoints->pushBack(newPoint[2]);
return index;
}
void
unsigned int
insertPointWithCheck(const std::vector<double> & newPoint,
std::map<std::vector<double>, unsigned int, PointComparison> & coordToIdMap,
const shared_ptr<XdmfTopology> & newConnectivity,
......@@ -108,10 +106,12 @@ namespace {
if(iter == coordToIdMap.end()) {
// Not inserted before
coordToIdMap[newPoint] = newPoints->getSize() / 3;;
insertPointWithoutCheck(newPoint, newConnectivity, newPoints);
return insertPointWithoutCheck(newPoint, newConnectivity, newPoints);
}
else {
newConnectivity->pushBack(iter->second);
const unsigned int index = iter->second;
newConnectivity->pushBack(index);
return index;
}
}
......@@ -127,6 +127,7 @@ namespace {
shared_ptr<XdmfUnstructuredGrid>
convert(const shared_ptr<XdmfUnstructuredGrid> gridToConvert,
const shared_ptr<const XdmfTopologyType> topologyType,
const shared_ptr<XdmfHeavyDataWriter> heavyDataWriter) const
{
shared_ptr<XdmfUnstructuredGrid> toReturn =
......@@ -223,1576 +224,494 @@ namespace {
};
class HexahedronToHexahedron_64 : public Converter {
template <unsigned int ORDER, bool ISSPECTRAL>
class HexahedronToHighOrderHexahedron : public Converter {
public:
HexahedronToHexahedron_64()
{
}
virtual ~HexahedronToHexahedron_64()
HexahedronToHighOrderHexahedron()
{
}
virtual void
computeInteriorPoints(std::vector<double> & leftPoint,
std::vector<double> & rightPoint,
const std::vector<double> & point1,
const std::vector<double> & point2) const
virtual ~HexahedronToHighOrderHexahedron()
{
this->computeLeftPoint(leftPoint, point1, point2);
this->computeRightPoint(rightPoint, point1, point2);
}
virtual void
computeLeftPoint(std::vector<double> & leftPoint,
const std::vector<double> & point1,
const std::vector<double> & point2) const
void
calculateIntermediatePoint(std::vector<double> & result,
const std::vector<double> & point1,
const std::vector<double> & point2,
int index,
bool spectral) const
{
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]);
const double scalar = points[index];
for (int i=0; i<3; i++)
result[i] = point1[i]+scalar*(point2[i]-point1[i]);
}
virtual void
computeRightPoint(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]);
}
shared_ptr<XdmfUnstructuredGrid>
convert(const shared_ptr<XdmfUnstructuredGrid> gridToConvert,
const shared_ptr<const XdmfTopologyType> topologyType,
const shared_ptr<XdmfHeavyDataWriter> heavyDataWriter) const
{
shared_ptr<XdmfUnstructuredGrid> toReturn =
XdmfUnstructuredGrid::New();
shared_ptr<XdmfUnstructuredGrid> toReturn = XdmfUnstructuredGrid::New();
toReturn->setName(gridToConvert->getName());
shared_ptr<XdmfGeometry> toReturnGeometry =
toReturn->getGeometry();
toReturnGeometry->setType(gridToConvert->getGeometry()->getType());
toReturnGeometry->initialize(gridToConvert->getGeometry()->getArrayType(),
gridToConvert->getGeometry()->getSize());
shared_ptr<XdmfGeometry> geometry = gridToConvert->getGeometry();
shared_ptr<XdmfGeometry> toReturnGeometry = toReturn->getGeometry();
toReturnGeometry->setType(geometry->getType());
toReturnGeometry->initialize(geometry->getArrayType());
bool releaseGeometry = false;
if(!gridToConvert->getGeometry()->isInitialized()) {
gridToConvert->getGeometry()->read();
if(!geometry->isInitialized()) {
geometry->read();
releaseGeometry = true;
}
// Copy all geometry values from old grid into new grid because we are
// keeping all old points.
toReturnGeometry->insert(0,
gridToConvert->getGeometry(),
0,
gridToConvert->getGeometry()->getSize());
if(releaseGeometry) {
gridToConvert->getGeometry()->release();
}
shared_ptr<XdmfTopology> topology = gridToConvert->getTopology();
shared_ptr<XdmfTopology> toReturnTopology = toReturn->getTopology();
shared_ptr<XdmfTopology> toReturnTopology =
toReturn->getTopology();
toReturnTopology->setType(XdmfTopologyType::Hexahedron_64());
toReturnTopology->initialize(gridToConvert->getTopology()->getArrayType());
toReturnTopology->reserve(64 * gridToConvert->getTopology()->getNumberElements());
toReturn->getTopology()->setType(topologyType);
toReturnTopology->initialize(topology->getArrayType());
toReturnTopology->reserve(mNodesPerElement *
topology->getNumberElements());
bool releaseTopology = false;
if(!gridToConvert->getTopology()->isInitialized()) {
gridToConvert->getTopology()->read();
if(!topology->isInitialized()) {
topology->read();
releaseTopology = true;
}
std::vector<double> leftPoint(3);
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.
for(int j=0; j<8; ++j) {
toReturnGeometry->getValues(gridToConvert->getTopology()->getValue<unsigned int>(8*i + j) * 3,
&localNodes[j][0], 3);
std::map<unsigned int, unsigned int> 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<double> elementCorner0(3);
std::vector<double> elementCorner1(3);
std::vector<double> elementCorner2(3);
std::vector<double> elementCorner3(3);
std::vector<double> elementCorner4(3);
std::vector<double> elementCorner5(3);
std::vector<double> elementCorner6(3);
std::vector<double> elementCorner7(3);
std::vector<double> planeCorner0(3);
std::vector<double> planeCorner1(3);
std::vector<double> planeCorner2(3);
std::vector<double> planeCorner3(3);
std::vector<double> lineEndPoint0(3);
std::vector<double> lineEndPoint1(3);
std::vector<double> point(3);
unsigned int offset = 0;
for(unsigned int elem = 0; elem<topology->getNumberElements(); ++elem) {
//
// get indices of coner vertices of the element
//
zeroIndex = topology->getValue<unsigned int>(offset++);
oneIndex = topology->getValue<unsigned int>(offset++);
twoIndex = topology->getValue<unsigned int>(offset++);
threeIndex = topology->getValue<unsigned int>(offset++);
fourIndex = topology->getValue<unsigned int>(offset++);
fiveIndex = topology->getValue<unsigned int>(offset++);
sixIndex = topology->getValue<unsigned int>(offset++);
sevenIndex = topology->getValue<unsigned int>(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);
// loop over i, j, k directions of element isolation i, j, and
// k planes
for(unsigned int i=0; i<mNodesPerEdge; ++i){
// calculate corners of i plane
calculateIntermediatePoint(planeCorner0,
elementCorner0,
elementCorner1,
i,
true);
calculateIntermediatePoint(planeCorner1,
elementCorner4,
elementCorner5,
i,
true);
calculateIntermediatePoint(planeCorner2,
elementCorner3,
elementCorner2,
i,
true);
calculateIntermediatePoint(planeCorner3,
elementCorner7,
elementCorner6,
i,
true);
for(unsigned int j=0; j<mNodesPerEdge; ++j) {
// calculate endpoints of j slice of i plane
calculateIntermediatePoint(lineEndPoint0,
planeCorner0,
planeCorner2,
j,
true);
calculateIntermediatePoint(lineEndPoint1,
planeCorner1,
planeCorner3,
j,
true);
for(unsigned int k=0; k<mNodesPerEdge; ++k) {
// calculate point to add to mesh
calculateIntermediatePoint(point,
lineEndPoint0,
lineEndPoint1,
k,
true);
if((i == 0 || i == ORDER) ||
(j == 0 || j == ORDER) ||
(k == 0 || k == ORDER)) {
unsigned int newIndex =
this->insertPointWithCheck(point,
coordToIdMap,
toReturnTopology,
toReturnGeometry);
if((i == 0 || i == ORDER) &&
(j == 0 || j == ORDER) &&
(k == 0 || k == ORDER)) {
if(i == 0) {
if(j == 0) {
if(k == 0) {
oldIdToNewId[zeroIndex] = newIndex;
}
else {
oldIdToNewId[fourIndex] = newIndex;
}
}
else if(k == 0) {
oldIdToNewId[threeIndex] = newIndex;
}
else {
oldIdToNewId[sevenIndex] = newIndex;
}
}
else {
if(j == 0) {
if(k == 0) {
oldIdToNewId[oneIndex] = newIndex;
}
else {
oldIdToNewId[fiveIndex] = newIndex;
}
}
else if(k == 0) {
oldIdToNewId[twoIndex] = newIndex;
}
else {
oldIdToNewId[sixIndex] = newIndex;
}
}
}
}
else {
this->insertPointWithoutCheck(point,
toReturnTopology,
toReturnGeometry);
}
}
}
}
// Add old connectivity information to newConnectivity.
toReturnTopology->resize(toReturnTopology->getSize() + 8, 0);
toReturnTopology->insert(64*i, gridToConvert->getTopology(), 8*i, 8);