Commit 453fc041 authored by Dave Moxey's avatar Dave Moxey

Start output of high-order boundary surfaces

parent 41670258
......@@ -457,7 +457,8 @@ public:
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType edgeType,
int coordDim,
int &id)
int &id,
bool justConfig = false)
{
ASSERTL0(false,
"This function should be implemented at a shape level.");
......
......@@ -249,6 +249,7 @@ public:
const int nTriPts = nPoints * (nPoints + 1) / 2;
const int nTriIntPts = (nPoints - 3) * (nPoints - 2) / 2;
m_faceNodes.resize(nTriIntPts);
std::cout << "num tri int pts: " << nTriIntPts << " " << nTriPts << " " << nPoints << std::endl;
for (int i = 3 + 3*(nPoints-2), cnt = 0; i < nTriPts; ++i, ++cnt)
{
......
......@@ -202,8 +202,10 @@ void Hexahedron::MakeOrder(int order,
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType pType,
int coordDim,
int &id)
int &id,
bool justConfig)
{
m_curveType = pType;
m_conf.m_order = order;
m_volumeNodes.clear();
......@@ -215,7 +217,11 @@ void Hexahedron::MakeOrder(int order,
m_conf.m_faceNodes = true;
m_conf.m_volumeNodes = true;
m_curveType = pType;
if (justConfig)
{
return;
}
int nPoints = order + 1;
StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
......
......@@ -84,7 +84,8 @@ public:
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType pType,
int coordDim,
int &id);
int &id,
bool justConfig = false);
NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf);
NEKMESHUTILS_EXPORT virtual int GetFaceVertex(int i, int j)
......
......@@ -109,6 +109,67 @@ SpatialDomains::GeometrySharedPtr Line::GetGeom(int coordDim)
return ret;
}
void Line::MakeOrder(int order,
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType pType,
int coordDim,
int &id,
bool justConfig)
{
m_conf.m_order = order;
m_curveType = pType;
m_conf.m_volumeNodes = false;
m_volumeNodes.clear();
// Lines of order == 1 have no interior volume points.
if (order == 1)
{
m_conf.m_faceNodes = false;
return;
}
m_conf.m_faceNodes = true;
if (justConfig)
{
return;
}
int nPoints = order + 1;
StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
Array<OneD, NekDouble> px;
LibUtilities::PointsKey pKey(nPoints, pType);
ASSERTL1(pKey.GetPointsDim() == 1, "Points distribution must be 1D");
LibUtilities::PointsManager()[pKey]->GetPoints(px);
Array<OneD, Array<OneD, NekDouble> > phys(coordDim);
for (int i = 0; i < coordDim; ++i)
{
phys[i] = Array<OneD, NekDouble>(xmap->GetTotPoints());
xmap->BwdTrans(geom->GetCoeffs(i), phys[i]);
}
int nQuadIntPts = (nPoints - 2) * (nPoints - 2);
m_volumeNodes.resize(nQuadIntPts);
for (int i = 1, cnt = 0; i < nPoints-1; ++i)
{
Array<OneD, NekDouble> xp(1);
xp[0] = px[i];
Array<OneD, NekDouble> x(3, 0.0);
for (int k = 0; k < coordDim; ++k)
{
x[k] = xmap->PhysEvaluate(xp, phys[k]);
}
m_volumeNodes[cnt] = boost::shared_ptr<Node>(
new Node(id++, x[0], x[1], x[2]));
}
}
/**
* @brief Return the number of nodes defining a line.
*/
......
......@@ -68,7 +68,13 @@ public:
NEKMESHUTILS_EXPORT virtual SpatialDomains::GeometrySharedPtr GetGeom(
int coordDim);
NEKMESHUTILS_EXPORT virtual void MakeOrder(
int order,
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType pType,
int coordDim,
int &id,
bool justConfig = false);
NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf);
};
}
......
......@@ -96,27 +96,33 @@ void Mesh::MakeOrder(int order,
boost::unordered_map<int, SpatialDomains::Geometry2DSharedPtr> faceGeoms;
boost::unordered_map<int, SpatialDomains::GeometrySharedPtr> volGeoms;
// Decide on distribution of points to use for each shape type.
// Decide on distribution of points to use for each shape type based on the
// input we've been supplied.
std::map<LibUtilities::ShapeType, LibUtilities::PointsType> pTypes;
if (distType == LibUtilities::ePolyEvenlySpaced)
{
pTypes[LibUtilities::eSegment] = LibUtilities::ePolyEvenlySpaced;
pTypes[LibUtilities::eTriangle] = LibUtilities::eNodalTriEvenlySpaced;
pTypes[LibUtilities::eSegment] = LibUtilities::ePolyEvenlySpaced;
pTypes[LibUtilities::eTriangle] = LibUtilities::eNodalTriEvenlySpaced;
pTypes[LibUtilities::eQuadrilateral] = LibUtilities::ePolyEvenlySpaced;
pTypes[LibUtilities::eTetrahedron] = LibUtilities::eNodalTetEvenlySpaced;
pTypes[LibUtilities::ePrism] = LibUtilities::eNodalPrismEvenlySpaced;
pTypes[LibUtilities::eHexahedron] = LibUtilities::ePolyEvenlySpaced;
pTypes[LibUtilities::eTetrahedron] =
LibUtilities::eNodalTetEvenlySpaced;
pTypes[LibUtilities::ePrism] = LibUtilities::eNodalPrismEvenlySpaced;
pTypes[LibUtilities::eHexahedron] = LibUtilities::ePolyEvenlySpaced;
}
else if (distType == LibUtilities::eGaussLobattoLegendre)
{
pTypes[LibUtilities::eSegment] = LibUtilities::eGaussLobattoLegendre;
pTypes[LibUtilities::eTriangle] = LibUtilities::eNodalTriElec;
pTypes[LibUtilities::eQuadrilateral] = LibUtilities::eGaussLobattoLegendre;
pTypes[LibUtilities::eTetrahedron] = LibUtilities::eNodalTetElec;
// Prism still to do.
pTypes[LibUtilities::eHexahedron] = LibUtilities::eGaussLobattoLegendre;
pTypes[LibUtilities::eSegment] = LibUtilities::eGaussLobattoLegendre;
pTypes[LibUtilities::eTriangle] = LibUtilities::eNodalTriElec;
pTypes[LibUtilities::eQuadrilateral] =
LibUtilities::eGaussLobattoLegendre;
pTypes[LibUtilities::eTetrahedron] = LibUtilities::eNodalTetElec;
pTypes[LibUtilities::eHexahedron] = LibUtilities::eGaussLobattoLegendre;
}
// Begin by generating Nektar++ geometry objects for edges, faces and
// elements so that we don't affect any neighbouring elements in the mesh as
// we process each element.
for(eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++)
{
SpatialDomains::Geometry1DSharedPtr geom =
......@@ -144,6 +150,8 @@ void Mesh::MakeOrder(int order,
boost::unordered_set<int> processedEdges, processedFaces, processedVolumes;
// Call MakeOrder with our generated geometries on each edge to fill in edge
// interior nodes.
for(eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++)
{
int edgeId = (*eit)->m_id;
......@@ -158,6 +166,8 @@ void Mesh::MakeOrder(int order,
processedEdges.insert(edgeId);
}
// Call MakeOrder with our generated geometries on each face to fill in face
// interior nodes.
for(fit = m_faceSet.begin(); fit != m_faceSet.end(); fit++)
{
int faceId = (*fit)->m_id;
......@@ -174,6 +184,40 @@ void Mesh::MakeOrder(int order,
processedFaces.insert(faceId);
}
// Copy curvature into boundary conditions
for (int i = 0; i < m_element[1].size(); ++i)
{
ElementSharedPtr el = m_element[1][i];
EdgeSharedPtr edge = el->GetEdgeLink();
if (!edge)
{
continue;
}
// Copy face curvature
el->SetVolumeNodes(edge->m_edgeNodes);
el->MakeOrder(order, SpatialDomains::GeometrySharedPtr(),
pTypes[el->GetConf().m_e], m_spaceDim, id, true);
}
for (int i = 0; i < m_element[2].size(); ++i)
{
ElementSharedPtr el = m_element[2][i];
FaceSharedPtr face = el->GetFaceLink();
if (!face)
{
continue;
}
// Copy face curvature
el->MakeOrder(order, SpatialDomains::GeometrySharedPtr(),
pTypes[el->GetConf().m_e], m_spaceDim, id, true);
el->SetVolumeNodes(face->m_faceNodes);
}
// Finally, fill in volumes.
const int nElmt = m_element[m_expDim].size();
for (int i = 0; i < nElmt; ++i)
{
......
......@@ -281,9 +281,11 @@ void Prism::MakeOrder(int order,
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType pType,
int coordDim,
int &id)
int &id,
bool justConfig)
{
m_conf.m_order = order;
m_curveType = pType;
m_volumeNodes.clear();
if (order == 1)
......@@ -300,7 +302,11 @@ void Prism::MakeOrder(int order,
m_conf.m_faceNodes = true;
m_conf.m_volumeNodes = true;
m_curveType = pType;
if (justConfig)
{
return;
}
int nPoints = order + 1;
StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
......
......@@ -84,7 +84,9 @@ public:
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType pType,
int coordDim,
int &id);
int &id,
bool justConfig = false);
NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf);
NEKMESHUTILS_EXPORT virtual int GetFaceVertex(int i, int j)
{
......
......@@ -138,24 +138,25 @@ void Quadrilateral::MakeOrder(int order,
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType pType,
int coordDim,
int &id)
int &id,
bool justConfig)
{
m_conf.m_order = order;
m_conf.m_order = order;
m_curveType = pType;
m_conf.m_volumeNodes = false;
m_volumeNodes.clear();
// Quadrilaterals of order == 1 have no interior volume points.
if (order == 1)
{
m_conf.m_volumeNodes = m_conf.m_faceNodes = false;
m_conf.m_faceNodes = false;
return;
}
m_conf.m_faceNodes = true;
m_conf.m_volumeNodes = false;
m_curveType = pType;
m_conf.m_faceNodes = true;
// Quadrilaterals of order < 2 have no interior volume points.
if (order < 2)
if (justConfig)
{
m_volumeNodes.clear();
return;
}
......
......@@ -84,7 +84,8 @@ public:
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType pType,
int coordDim,
int &id);
int &id,
bool justConfig = false);
NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf);
};
......
......@@ -257,9 +257,11 @@ void Tetrahedron::MakeOrder(int order,
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType pType,
int coordDim,
int &id)
int &id,
bool justConfig)
{
m_conf.m_order = order;
m_curveType = pType;
m_volumeNodes.clear();
if (order == 1 || order == 2)
......@@ -279,10 +281,14 @@ void Tetrahedron::MakeOrder(int order,
return;
}
m_curveType = pType;
m_conf.m_faceNodes = true;
m_conf.m_volumeNodes = true;
if (justConfig)
{
return;
}
int nPoints = order + 1;
StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
......
......@@ -79,7 +79,8 @@ public:
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType pType,
int coordDim,
int &id);
int &id,
bool justConfig = false);
NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf);
NEKMESHUTILS_EXPORT virtual int GetFaceVertex(int i, int j)
......
......@@ -173,12 +173,25 @@ void Triangle::MakeOrder(int order,
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType pType,
int coordDim,
int &id)
int &id,
bool justConfig)
{
m_conf.m_order = order;
m_curveType = pType;
m_conf.m_volumeNodes = false;
m_volumeNodes.clear();
// Triangles of order < 3 have no interior volume points.
if (order < 3)
if (order == 1 || order == 2)
{
m_conf.m_faceNodes = false;
return;
}
m_conf.m_faceNodes = true;
if (justConfig)
{
m_volumeNodes.clear();
return;
}
......@@ -218,7 +231,6 @@ void Triangle::MakeOrder(int order,
new Node(id++, x[0], x[1], x[2]));
}
m_curveType = pType;
m_conf.m_order = order;
m_conf.m_faceNodes = true;
m_conf.m_volumeNodes = false;
......
......@@ -214,7 +214,8 @@ public:
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType pType,
int coordDim,
int &id);
int &id,
bool justConfig = false);
NEKMESHUTILS_EXPORT static unsigned int GetNumNodes(ElmtConfig pConf);
};
......
......@@ -1679,15 +1679,15 @@ std::map<unsigned int, ElmtConfig> InputGmsh::GenElmMap()
tmp[ 4] = ElmtConfig(eTetrahedron, 1, false, false);
tmp[ 5] = ElmtConfig(eHexahedron, 1, false, false);
tmp[ 6] = ElmtConfig(ePrism, 1, false, false);
tmp[ 7] = ElmtConfig(ePyramid, 1, true, true);
tmp[ 8] = ElmtConfig(eSegment, 2, true, true);
tmp[ 9] = ElmtConfig(eTriangle, 2, true, true);
tmp[ 10] = ElmtConfig(eQuadrilateral, 2, true, true);
tmp[ 11] = ElmtConfig(eTetrahedron, 2, true, true);
tmp[ 7] = ElmtConfig(ePyramid, 1, false, false);
tmp[ 8] = ElmtConfig(eSegment, 2, true, false);
tmp[ 9] = ElmtConfig(eTriangle, 2, false, false);
tmp[ 10] = ElmtConfig(eQuadrilateral, 2, true, false);
tmp[ 11] = ElmtConfig(eTetrahedron, 2, false, false);
tmp[ 12] = ElmtConfig(eHexahedron, 2, true, true);
tmp[ 13] = ElmtConfig(ePrism, 2, true, true);
tmp[ 14] = ElmtConfig(ePyramid, 2, true, true);
tmp[ 15] = ElmtConfig(ePoint, 1, true, false);
tmp[ 13] = ElmtConfig(ePrism, 2, true, false);
tmp[ 14] = ElmtConfig(ePyramid, 2, true, false);
tmp[ 15] = ElmtConfig(ePoint, 1, false, false);
tmp[ 16] = ElmtConfig(eQuadrilateral, 2, false, false);
tmp[ 17] = ElmtConfig(eHexahedron, 2, false, false);
tmp[ 18] = ElmtConfig(ePrism, 2, false, false);
......
......@@ -212,11 +212,18 @@ void OutputGmsh::Process()
tags.push_back(nodeList[j]->m_id);
}
cout << "GOT " << nodeList.size() << " VERTS" << endl;
cout << "GOT " << edgeList.size() << " EDGES" << endl;
cout << "GOT " << faceList.size() << " FACES" << endl;
cout << "GOT " << volList.size() << " VOLS" << endl;
// Process edge-interior points
for (int j = 0; j < edgeList.size(); ++j)
{
nodeList = edgeList[j]->m_edgeNodes;
cout << "edge " << j << ": " << edgeList[j]->m_edgeNodes.size() << endl;
if (e->GetEdgeOrient(j, edgeList[j]) == StdRegions::eForwards)
{
for (int k = 0; k < nodeList.size(); ++k)
......@@ -269,6 +276,7 @@ void OutputGmsh::Process()
}
// Process volume nodes
cout << volList.size() << endl;
for (int j = 0; j < volList.size(); ++j)
{
tags.push_back(volList[j]->m_id);
......@@ -278,6 +286,12 @@ void OutputGmsh::Process()
vector<int> reordering = InputGmsh::CreateReordering(elmtType);
vector<int> inv(tags.size());
cout << "TAGS LENGTH: " << tags.size() << endl;
cout << "REORDERING LENGTH: " << reordering.size() << endl;
ASSERTL1(tags.size() == reordering.size(),
"Reordering map size not equal to element tags.");
for (int j = 0; j < tags.size(); ++j)
{
inv[reordering[j]] = j;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment