Commit 6a5d3eb6 authored by Michael Turner's avatar Michael Turner

add cad logic to make order

parent a9031e74
......@@ -95,6 +95,29 @@ NekDouble CADCurve::Length(NekDouble ti, NekDouble tf)
return System.Mass() / 1000.0;
}
void CADCurve::Loct(Array<OneD, NekDouble> &xyz, NekDouble &t)
{
Array<OneD, NekDouble> b = Bounds();
Handle(Geom_Curve) NewCurve = BRep_Tool::Curve(m_occEdge, b[0], b[1]);
gp_Pnt loc(xyz[0] * 1000.0, xyz[1] * 1000.0, xyz[2] * 1000.0);
GeomAPI_ProjectPointOnCurve locator(loc,NewCurve);
if (locator.NbPoints() == 0)
{
ASSERTL0(false,"failed");
}
else
{
t = locator.Parameter(0);
gp_Pnt tst = locator.NearestPoint();
NekDouble dis = tst.Distance(loc);
if(dis > 1e-6)
{
cout << "large curve projection: " << dis << endl;
}
}
}
Array<OneD, NekDouble> CADCurve::P(NekDouble t)
{
Array<OneD, NekDouble> location(3);
......
......@@ -154,6 +154,8 @@ public:
return m_mainVerts;
}
void Loct(Array<OneD, NekDouble> &xyz, NekDouble &t);
private:
/// OpenCascade object of the curve.
BRepAdaptor_Curve m_occCurve;
......
......@@ -48,6 +48,7 @@
#include <BRepAdaptor_Curve.hxx>
#include <BRepAdaptor_Surface.hxx>
#include <GeomAPI_ProjectPointOnSurf.hxx>
#include <GeomAPI_ProjectPointOnCurve.hxx>
#include <GeomAbs_SurfaceType.hxx>
#include <BRepTools.hxx>
#include <BRep_Tool.hxx>
......
SET(NEKMESHUTILS_SOURCES
MeshElements/Mesh.cpp
MeshElements/Element.cpp
MeshElements/Edge.cpp
MeshElements/Face.cpp
MeshElements/BooleanOperators.cpp
MeshElements/Point.cpp
MeshElements/Line.cpp
......
......@@ -42,7 +42,6 @@
#include <SpatialDomains/SegGeom.h>
#include <NekMeshUtils/NekMeshUtilsDeclspec.h>
#include <NekMeshUtils/MeshElements/Element.h>
#include <NekMeshUtils/MeshElements/Node.h>
namespace Nektar
......@@ -72,6 +71,7 @@ public:
{
#ifdef NEKTAR_USE_MESHGEN
onCurve = false;
onSurf = false;
#endif
}
......@@ -81,6 +81,7 @@ public:
{
#ifdef NEKTAR_USE_MESHGEN
onCurve = false;
onSurf = false;
#endif
}
......@@ -114,103 +115,17 @@ public:
/// Creates a Nektar++ string listing the coordinates of all the
/// nodes.
NEKMESHUTILS_EXPORT std::string GetXmlCurveString() const
{
std::vector<NodeSharedPtr> nodeList;
GetCurvedNodes(nodeList);
std::stringstream s;
std::string str;
// put them into a string and return
for (int k = 0; k < nodeList.size(); ++k)
{
s << std::scientific << std::setprecision(8) << " "
<< nodeList[k]->m_x << " " << nodeList[k]->m_y << " "
<< nodeList[k]->m_z << " ";
}
return s.str();
}
NEKMESHUTILS_EXPORT std::string GetXmlCurveString();
/// Generate a SpatialDomains::SegGeom object for this edge.
NEKMESHUTILS_EXPORT SpatialDomains::SegGeomSharedPtr GetGeom(int coordDim)
{
static boost::mutex io_mutex;
boost::mutex::scoped_lock lock(io_mutex);
// Create edge vertices.
SpatialDomains::PointGeomSharedPtr p[2];
SpatialDomains::SegGeomSharedPtr ret;
p[0] = m_n1->GetGeom(coordDim);
p[1] = m_n2->GetGeom(coordDim);
// Create a curve if high-order information exists.
if (m_edgeNodes.size() > 0)
{
SpatialDomains::CurveSharedPtr c =
MemoryManager<SpatialDomains::Curve>::AllocateSharedPtr(
m_id, m_curveType);
c->m_points.push_back(p[0]);
for (int i = 0; i < m_edgeNodes.size(); ++i)
{
c->m_points.push_back(m_edgeNodes[i]->GetGeom(coordDim));
}
c->m_points.push_back(p[1]);
ret = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(
m_id, coordDim, p, c);
}
else
{
ret = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(
m_id, coordDim, p);
}
return ret;
}
NEKMESHUTILS_EXPORT SpatialDomains::SegGeomSharedPtr GetGeom(int coordDim);
/// Make this edge an order @p order edge. @see Element::MakeOrder.
void MakeOrder(int order,
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType edgeType,
int coordDim,
int &id)
{
int nPoints = order + 1;
StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
Array<OneD, NekDouble> edgePoints;
LibUtilities::PointsKey edgeKey(nPoints, edgeType);
LibUtilities::PointsManager()[edgeKey]->GetPoints(edgePoints);
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]);
}
m_edgeNodes.resize(nPoints - 2);
for (int i = 1; i < nPoints - 1; ++i)
{
Array<OneD, NekDouble> x(3, 0.0);
for (int j = 0; j < coordDim; ++j)
{
x[j] = xmap->PhysEvaluate(edgePoints + i, phys[j]);
}
m_edgeNodes[i-1] = boost::shared_ptr<Node>(
new Node(id++, x[0], x[1], x[2]));
}
m_curveType = edgeType;
}
int &id);
/// ID of edge.
unsigned int m_id;
......@@ -230,6 +145,9 @@ public:
/// id of cad curve which edge lies on
int CADCurveId;
CADCurveSharedPtr CADCurve;
bool onSurf;
int CADSurfId;
CADSurfSharedPtr CADSurf;
#endif
private:
......
......@@ -41,6 +41,7 @@
#include <NekMeshUtils/NekMeshUtilsDeclspec.h>
#include <NekMeshUtils/MeshElements/Edge.h>
#include <NekMeshUtils/MeshElements/Node.h>
namespace Nektar
{
......@@ -112,361 +113,22 @@ public:
/// Assemble a list of nodes on curved face
NEKMESHUTILS_EXPORT void GetCurvedNodes(
std::vector<NodeSharedPtr> &nodeList) const
{
// Treat 2D point distributions differently to 3D.
if (m_curveType == LibUtilities::eNodalTriFekete ||
m_curveType == LibUtilities::eNodalTriEvenlySpaced ||
m_curveType == LibUtilities::eNodalTriElec)
{
int n = m_edgeList[0]->GetNodeCount();
nodeList.insert(
nodeList.end(), m_vertexList.begin(), m_vertexList.end());
for (int k = 0; k < m_edgeList.size(); ++k)
{
nodeList.insert(nodeList.end(),
m_edgeList[k]->m_edgeNodes.begin(),
m_edgeList[k]->m_edgeNodes.end());
if (m_edgeList[k]->m_n1 != m_vertexList[k])
{
// If edge orientation is reversed relative to node
// ordering, we need to reverse order of nodes.
std::reverse(nodeList.begin() + 3 + k * (n - 2),
nodeList.begin() + 3 + (k + 1) * (n - 2));
}
}
nodeList.insert(
nodeList.end(), m_faceNodes.begin(), m_faceNodes.end());
}
else
{
// Write out in 2D tensor product order.
ASSERTL0(m_vertexList.size() == 4,
"Face nodes of tensor product only supported "
"for quadrilaterals.");
int n = (int)sqrt((NekDouble)GetNodeCount());
nodeList.resize(n * n);
ASSERTL0(n * n == GetNodeCount(), "Wrong number of modes?");
// Write vertices
nodeList[0] = m_vertexList[0];
nodeList[n - 1] = m_vertexList[1];
nodeList[n * n - 1] = m_vertexList[2];
nodeList[n * (n - 1)] = m_vertexList[3];
// Write edge-interior
int skips[4][2] = {
{0, 1}, {n - 1, n}, {n * n - 1, -1}, {n * (n - 1), -n}};
for (int i = 0; i < 4; ++i)
{
bool reverseEdge = m_edgeList[i]->m_n1 == m_vertexList[i];
if (!reverseEdge)
{
for (int j = 1; j < n - 1; ++j)
{
nodeList[skips[i][0] + j * skips[i][1]] =
m_edgeList[i]->m_edgeNodes[n - 2 - j];
}
}
else
{
for (int j = 1; j < n - 1; ++j)
{
nodeList[skips[i][0] + j * skips[i][1]] =
m_edgeList[i]->m_edgeNodes[j - 1];
}
}
}
// Write interior
for (int i = 1; i < n - 1; ++i)
{
for (int j = 1; j < n - 1; ++j)
{
nodeList[i * n + j] =
m_faceNodes[(i - 1) * (n - 2) + (j - 1)];
}
}
}
}
std::vector<NodeSharedPtr> &nodeList);
/// Generates a string listing the coordinates of all nodes
/// associated with this face.
NEKMESHUTILS_EXPORT std::string GetXmlCurveString() const
{
std::stringstream s;
std::string str;
std::vector<NodeSharedPtr> nodeList;
// assemble listof nodes
GetCurvedNodes(nodeList);
// put them into a string
for (int k = 0; k < nodeList.size(); ++k)
{
s << std::scientific << std::setprecision(8) << " "
<< nodeList[k]->m_x << " " << nodeList[k]->m_y << " "
<< nodeList[k]->m_z << " ";
}
return s.str();
}
NEKMESHUTILS_EXPORT std::string GetXmlCurveString();
/// Make this face an order @p order face. @see Element::MakeOrder.
void MakeOrder(int order,
SpatialDomains::GeometrySharedPtr geom,
LibUtilities::PointsType pType,
int coordDim,
int &id)
{
if (m_vertexList.size() == 3)
{
// Triangles of order < 3 have no interior volume points.
if (order < 3)
{
m_faceNodes.clear();
return;
}
int nPoints = order + 1;
StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
Array<OneD, NekDouble> px, py;
LibUtilities::PointsKey pKey(nPoints, pType);
ASSERTL1(pKey.GetPointsDim() == 2, "Points distribution must be 2D");
LibUtilities::PointsManager()[pKey]->GetPoints(px, py);
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]);
}
const int nTriPts = nPoints * (nPoints + 1) / 2;
const int nTriIntPts = (nPoints - 3) * (nPoints - 2) / 2;
m_faceNodes.resize(nTriIntPts);
for (int i = 3 + 3*(nPoints-2), cnt = 0; i < nTriPts; ++i, ++cnt)
{
Array<OneD, NekDouble> xp(2);
xp[0] = px[i];
xp[1] = py[i];
Array<OneD, NekDouble> x(3, 0.0);
for (int j = 0; j < coordDim; ++j)
{
x[j] = xmap->PhysEvaluate(xp, phys[j]);
}
m_faceNodes[cnt] = boost::shared_ptr<Node>(
new Node(id++, x[0], x[1], x[2]));
}
m_curveType = pType;
}
else if (m_vertexList.size() == 4)
{
// Quads of order < 2 have no interior volume points.
if (order < 2)
{
m_faceNodes.clear();
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_faceNodes.resize(nQuadIntPts);
for (int i = 1, cnt = 0; i < nPoints-1; ++i)
{
for (int j = 1; j < nPoints-1; ++j, ++cnt)
{
Array<OneD, NekDouble> xp(2);
xp[0] = px[j];
xp[1] = px[i];
Array<OneD, NekDouble> x(3, 0.0);
for (int k = 0; k < coordDim; ++k)
{
x[k] = xmap->PhysEvaluate(xp, phys[k]);
}
m_faceNodes[cnt] = boost::shared_ptr<Node>(
new Node(id++, x[0], x[1], x[2]));
}
}
m_curveType = pType;
}
else
{
ASSERTL0(false, "Unknown number of vertices");
}
}
int &id);
/// Generate either SpatialDomains::TriGeom or
/// SpatialDomains::QuadGeom for this element.
NEKMESHUTILS_EXPORT SpatialDomains::Geometry2DSharedPtr GetGeom(int coordDim)
{
int nEdge = m_edgeList.size();
SpatialDomains::SegGeomSharedPtr edges[4];
SpatialDomains::Geometry2DSharedPtr ret;
StdRegions::Orientation edgeo[4];
for (int i = 0; i < nEdge; ++i)
{
edges[i] = m_edgeList[i]->GetGeom(coordDim);
}
for (int i = 0; i < nEdge; ++i)
{
edgeo[i] = m_edgeList[i]->m_n1 == m_vertexList[i]
? StdRegions::eForwards
: StdRegions::eBackwards;
}
if (m_faceNodes.size() > 0)
{
if (nEdge == 3)
{
SpatialDomains::CurveSharedPtr c =
MemoryManager<SpatialDomains::Curve>::AllocateSharedPtr(
m_id, m_curveType);
for (int j = 0; j < m_vertexList.size(); j++)
{
c->m_points.push_back(m_vertexList[j]->GetGeom(coordDim));
}
for (int j = 0; j < nEdge; j++)
{
std::vector<NodeSharedPtr> ed = m_edgeList[j]->m_edgeNodes;
if (edgeo[j] == StdRegions::eBackwards)
{
for (int k = ed.size() - 1; k >= 0; k--)
{
c->m_points.push_back(ed[k]->GetGeom(coordDim));
}
}
else
{
for (int k = 0; k < ed.size(); k++)
{
c->m_points.push_back(ed[k]->GetGeom(coordDim));
}
}
}
for (int j = 0; j < m_faceNodes.size(); j++)
{
c->m_points.push_back(m_faceNodes[j]->GetGeom(coordDim));
}
ret = MemoryManager<SpatialDomains::TriGeom>::AllocateSharedPtr(
m_id, edges, edgeo, c);
}
else
{
SpatialDomains::CurveSharedPtr c =
MemoryManager<SpatialDomains::Curve>::AllocateSharedPtr(
m_id, m_curveType);
ASSERTL0(m_vertexList.size() == 4,
"Face nodes of tensor product only supported "
"for quadrilaterals.");
int n = (int)sqrt((NekDouble)GetNodeCount());
std::vector<NodeSharedPtr> tmp(n * n);
ASSERTL0(n * n == GetNodeCount(), "Wrong number of modes?");
// Write vertices
tmp[0] = m_vertexList[0];
tmp[n - 1] = m_vertexList[1];
tmp[n * n - 1] = m_vertexList[2];
tmp[n * (n - 1)] = m_vertexList[3];
// Write edge-interior
int skips[4][2] = {
{0, 1}, {n - 1, n}, {n * n - 1, -1}, {n * (n - 1), -n}};
for (int i = 0; i < 4; ++i)
{
bool reverseEdge = edgeo[i] == StdRegions::eBackwards;
if (reverseEdge)
{
for (int j = 1; j < n - 1; ++j)
{
tmp[skips[i][0] + j * skips[i][1]] =
m_edgeList[i]->m_edgeNodes[n - 2 - j];
}
}
else
{
for (int j = 1; j < n - 1; ++j)
{
tmp[skips[i][0] + j * skips[i][1]] =
m_edgeList[i]->m_edgeNodes[j - 1];
}
}
}
// Write interior
for (int i = 1; i < n - 1; ++i)
{
for (int j = 1; j < n - 1; ++j)
{
tmp[i * n + j] =
m_faceNodes[(i - 1) * (n - 2) + (j - 1)];
}
}
for (int k = 0; k < tmp.size(); ++k)
{
c->m_points.push_back(tmp[k]->GetGeom(coordDim));
}
ret =
MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(
m_id, edges, edgeo, c);
}
}
else
{
if (nEdge == 3)
{
ret = MemoryManager<SpatialDomains::TriGeom>::AllocateSharedPtr(
m_id, edges, edgeo);
}
else
{
ret =
MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(
m_id, edges, edgeo);
}
}
return ret;
}
NEKMESHUTILS_EXPORT SpatialDomains::Geometry2DSharedPtr GetGeom(int coordDim);
/// ID of the face.
unsigned int m_id;
......
......@@ -595,6 +595,7 @@ void InputNekpp::Process()
for(int j = 0; j < (*it)->m_edgeNodes.size(); j++)
{
esit = edgeToString.find(pair<int,int>((*it)->m_id,j));
int surf = 0;
if(esit != edgeToString.end())
{
ct++;
......@@ -609,6 +610,7 @@ void InputNekpp::Process()
else if(esit->second.second[i].type == "S")
{
int s = esit->second.second[i].id;
surf = s;
Array<OneD,NekDouble> uv(2);
uv[0] = esit->second.second[i].u;
uv[1] = esit->second.second[i].v;
......@@ -626,6 +628,12 @@ void InputNekpp::Process()
(*it)->CADCurveId = esit->second.first;
(*it)->CADCurve = m_mesh->m_cad->GetCurve(esit->second.first);
}
else if(surf > 0)