Commit c140cd2d authored by Michael Turner's avatar Michael Turner
Browse files

Merge branch 'feature/non-O-BL' into 'master'

Non-O BL meshing in 2D

See merge request !757
parents 47302e60 81a84a8a
......@@ -6,6 +6,7 @@ v4.5.0
**NekMesh**:
- Add periodic boundary condition meshing in 2D (!733)
- Adjust boundary layer thickness in corners in 2D (!739)
- Add non-O BL meshing in 2D (!757)
**Library**
- Added in sum factorisation version for pyramid expansions and orthogonal
......
......@@ -83,6 +83,12 @@ void Generator2D::Process()
ParseUtils::GenerateSeqVector(m_config["blcurves"].as<string>().c_str(),
m_blCurves);
// find the ends of the BL curves
if (m_config["blcurves"].beenSet)
{
FindBLEnds();
}
// linear mesh all curves
for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); i++)
{
......@@ -97,6 +103,34 @@ void Generator2D::Process()
{
m_curvemeshes[i] =
MemoryManager<CurveMesh>::AllocateSharedPtr(i, m_mesh);
// Fheck if this curve is at an end of the BL
// If so, define an offset for the second node, corresponding to the
// BL thickness
if (m_blends.count(i))
{
vector<CADVertSharedPtr> vertices =
m_mesh->m_cad->GetCurve(i)->GetVertex();
Array<OneD, NekDouble> loc;
NekDouble t;
// offset needed at first node (or both)
if (m_blends[i] == 0 || m_blends[i] == 2)
{
loc = vertices[0]->GetLoc();
t = m_thickness.Evaluate(m_thickness_ID, loc[0], loc[1],
loc[2], 0.0);
m_curvemeshes[i]->SetOffset(0, t);
}
// offset needed at second node (or both)
if (m_blends[i] == 1 || m_blends[i] == 2)
{
loc = vertices[1]->GetLoc();
t = m_thickness.Evaluate(m_thickness_ID, loc[0], loc[1],
loc[2], 0.0);
m_curvemeshes[i]->SetOffset(1, t);
}
}
}
else
{
......@@ -125,6 +159,30 @@ void Generator2D::Process()
{
MakeBL(i);
}
// If the BL doesn't form closed loops, we need to remove the outside
// nodes from the curve meshes
for (map<unsigned, unsigned>::iterator ic = m_blends.begin();
ic != m_blends.end(); ++ic)
{
vector<NodeSharedPtr> nodes =
m_curvemeshes[ic->first]->GetMeshPoints();
if (ic->second == 0 || ic->second == 2)
{
nodes.erase(nodes.begin());
}
if (ic->second == 1 || ic->second == 2)
{
nodes.erase(nodes.end() - 1);
}
// Rebuild the curvemesh without the first node, the last node or
// both
m_curvemeshes[ic->first] =
MemoryManager<CurveMesh>::AllocateSharedPtr(ic->first, m_mesh,
nodes);
}
}
if (m_mesh->m_verbose)
......@@ -175,6 +233,67 @@ void Generator2D::Process()
Report();
}
void Generator2D::FindBLEnds()
{
// Set of CAD vertices
// Vertices of each curve are added to the set if not found and removed from
// the set if found
// This leaves us with a set of vertices that are at the end of BL open
// loops
set<CADVertSharedPtr> cadverts;
for (int it = 0; it < m_blCurves.size(); ++it)
{
vector<CADVertSharedPtr> vertices =
m_mesh->m_cad->GetCurve(m_blCurves[it])->GetVertex();
for (int iv = 0; iv < vertices.size(); ++iv)
{
set<CADVertSharedPtr>::iterator is = cadverts.find(vertices[iv]);
if (is != cadverts.end())
{
cadverts.erase(is);
}
else
{
cadverts.insert(vertices[iv]);
}
}
}
// Build m_blends based on the previously constructed set of vertices
// m_blends is a map of curve number (the curves right outside the BL open
// loops) to the offset node number: 0, 1 or 2 (for both)
for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); ++i)
{
if (find(m_blCurves.begin(), m_blCurves.end(), i) != m_blCurves.end())
{
continue;
}
vector<CADVertSharedPtr> vertices =
m_mesh->m_cad->GetCurve(i)->GetVertex();
for (int j = 0; j < 2; ++j)
{
if (!cadverts.count(vertices[j]))
{
continue;
}
if (m_blends.count(i))
{
m_blends[i] = 2;
}
else
{
m_blends[i] = j;
}
}
}
}
void Generator2D::MakeBLPrep()
{
if (m_mesh->m_verbose)
......@@ -184,10 +303,10 @@ void Generator2D::MakeBLPrep()
// identify the nodes which will become the boundary layer.
for (vector<unsigned>::iterator it = m_blCurves.begin();
it != m_blCurves.end(); ++it)
for (int it = 0; it < m_blCurves.size(); ++it)
{
vector<EdgeSharedPtr> localedges = m_curvemeshes[*it]->GetMeshEdges();
vector<EdgeSharedPtr> localedges =
m_curvemeshes[m_blCurves[it]]->GetMeshEdges();
for (int i = 0; i < localedges.size(); i++)
{
m_nodesToEdge[localedges[i]->m_n1].push_back(localedges[i]);
......@@ -200,12 +319,12 @@ void Generator2D::MakeBL(int faceid)
{
map<int, Array<OneD, NekDouble> > edgeNormals;
int eid = 0;
for (vector<unsigned>::iterator it = m_blCurves.begin();
it != m_blCurves.end(); ++it)
for (int it = 0; it < m_blCurves.size(); ++it)
{
CADOrientation::Orientation edgeo =
m_mesh->m_cad->GetCurve(*it)->GetOrienationWRT(faceid);
vector<EdgeSharedPtr> es = m_curvemeshes[*it]->GetMeshEdges();
m_mesh->m_cad->GetCurve(m_blCurves[it])->GetOrienationWRT(faceid);
vector<EdgeSharedPtr> es =
m_curvemeshes[m_blCurves[it]]->GetMeshEdges();
// on each !!!EDGE!!! calculate a normal
// always to the left unless edgeo is 1
// normal must be done in the parametric space (and then projected back)
......@@ -255,9 +374,32 @@ void Generator2D::MakeBL(int faceid)
map<NodeSharedPtr, vector<EdgeSharedPtr> >::iterator it;
for (it = m_nodesToEdge.begin(); it != m_nodesToEdge.end(); it++)
{
ASSERTL0(it->second.size() == 1 || it->second.size() == 2,
"weirdness, most likely bl_surfs are incorrect");
// If node at the end of the BL open loop, the "normal node" isn't
// constructed by computing a normal but found on the adjacent curve
if (it->second.size() == 1)
{
vector<pair<int, CADCurveSharedPtr> > curves =
it->first->GetCADCurves();
vector<EdgeSharedPtr> edges =
m_curvemeshes[curves[0].first]->GetMeshEdges();
vector<EdgeSharedPtr>::iterator ie =
find(edges.begin(), edges.end(), it->second[0]);
int rightCurve =
(ie == edges.end()) ? curves[0].first : curves[1].first;
vector<NodeSharedPtr> nodes =
m_curvemeshes[rightCurve]->GetMeshPoints();
nodeNormals[it->first] =
(nodes[0] == it->first) ? nodes[1] : nodes[nodes.size() - 2];
continue;
}
Array<OneD, NekDouble> n(3, 0.0);
ASSERTL0(it->second.size() == 2,
"wierdness, most likely bl_surfs are incorrect");
Array<OneD, NekDouble> n1 = edgeNormals[it->second[0]->m_id];
Array<OneD, NekDouble> n2 = edgeNormals[it->second[1]->m_id];
n[0] = (n1[0] + n2[0]) / 2.0;
......@@ -288,19 +430,20 @@ void Generator2D::MakeBL(int faceid)
nodeNormals[it->first] = nn;
}
for (vector<unsigned>::iterator it = m_blCurves.begin();
it != m_blCurves.end(); ++it)
for (int it = 0; it < m_blCurves.size(); ++it)
{
CADOrientation::Orientation edgeo =
m_mesh->m_cad->GetCurve(*it)->GetOrienationWRT(faceid);
vector<NodeSharedPtr> ns = m_curvemeshes[*it]->GetMeshPoints();
m_mesh->m_cad->GetCurve(m_blCurves[it])->GetOrienationWRT(faceid);
vector<NodeSharedPtr> ns =
m_curvemeshes[m_blCurves[it]]->GetMeshPoints();
vector<NodeSharedPtr> newNs;
for (int i = 0; i < ns.size(); i++)
{
newNs.push_back(nodeNormals[ns[i]]);
}
m_curvemeshes[*it] =
MemoryManager<CurveMesh>::AllocateSharedPtr(*it, m_mesh, newNs);
m_curvemeshes[m_blCurves[it]] =
MemoryManager<CurveMesh>::AllocateSharedPtr(m_blCurves[it], m_mesh,
newNs);
if (edgeo == CADOrientation::eBackwards)
{
reverse(ns.begin(), ns.end());
......@@ -346,10 +489,10 @@ void Generator2D::PeriodicPrep()
boost::split(lines, s, boost::is_any_of(":"));
for (vector<string>::iterator il = lines.begin(); il != lines.end(); ++il)
for (int il = 0; il < lines.size(); ++il)
{
vector<string> tmp;
boost::split(tmp, *il, boost::is_any_of(","));
boost::split(tmp, lines[il], boost::is_any_of(","));
ASSERTL0(tmp.size() == 2, "periodic pairs ill-defined");
......
......@@ -67,6 +67,8 @@ public:
virtual void Process();
private:
void FindBLEnds();
void MakeBLPrep();
void PeriodicPrep();
......@@ -84,6 +86,8 @@ private:
std::map<unsigned, unsigned> m_periodicPairs;
std::vector<unsigned int> m_blCurves;
/// map of curves and Bl ends: 0, 1 or 2 (for both)
std::map<unsigned, unsigned> m_blends;
LibUtilities::AnalyticExpressionEvaluator m_thickness;
int m_thickness_ID;
std::map<NodeSharedPtr, std::vector<EdgeSharedPtr> > m_nodesToEdge;
......
......@@ -52,6 +52,16 @@ void CurveMesh::Mesh()
int(m_curvelength / m_mesh->m_octree->GetMinDelta()) + 10;
ds = m_curvelength / (m_numSamplePoints - 1);
// compute the offset due to adjacent BLs
NekDouble totalOffset = 0.0;
for (map<unsigned, NekDouble>::iterator ie = m_endoffset.begin();
ie != m_endoffset.end(); ++ie)
{
totalOffset += ie->second;
}
ASSERTL0(m_curvelength > totalOffset,
"Boundary layers too thick for adjacent curve");
GetSampleFunction();
Ae = 0.0;
......@@ -63,12 +73,22 @@ void CurveMesh::Mesh()
Ne = round(Ae);
if (Ne + 1 < 2)
if (Ne + 1 < 2 + m_endoffset.size())
{
meshsvalue.resize(2);
Ne = 1 + m_endoffset.size();
meshsvalue.resize(Ne + 1);
meshsvalue[0] = 0.0;
meshsvalue[1] = m_curvelength;
Ne = 1;
if (m_endoffset.count(0))
{
meshsvalue[1] = m_endoffset[0];
}
if (m_endoffset.count(1))
{
meshsvalue[Ne - 1] = m_curvelength - m_endoffset[1];
}
}
else
{
......@@ -79,7 +99,19 @@ void CurveMesh::Mesh()
meshsvalue[0] = 0.0;
meshsvalue[Ne] = m_curvelength;
for (int i = 1; i < Ne; i++)
// force the second and/or the second to last point(s) if an offset is
// defined
if (m_endoffset.count(0))
{
meshsvalue[1] = m_endoffset[0];
}
if (m_endoffset.count(1))
{
meshsvalue[Ne - 1] = m_curvelength - m_endoffset[1];
}
for (int i = 1 + m_endoffset.count(0); i < Ne - m_endoffset.count(1);
i++)
{
int iterationcounter = 0;
bool iterate = true;
......@@ -108,7 +140,8 @@ void CurveMesh::Mesh()
Array<OneD, NekDouble> loc;
vector<CADVertSharedPtr> verts = m_cadcurve->GetVertex();
vector<pair<CADSurfSharedPtr, CADOrientation::Orientation> > s = m_cadcurve->GetAdjSurf();
vector<pair<CADSurfSharedPtr, CADOrientation::Orientation> > s =
m_cadcurve->GetAdjSurf();
NodeSharedPtr n = verts[0]->GetNode();
t = m_bounds[0];
......@@ -214,7 +247,7 @@ NekDouble CurveMesh::EvaluateDS(NekDouble s)
int a = 0;
int b = 0;
ASSERTL1(!(s < 0) && !(s > m_curvelength),"s out of bounds");
ASSERTL1(!(s < 0) && !(s > m_curvelength), "s out of bounds");
if (s == 0)
{
......@@ -253,7 +286,7 @@ NekDouble CurveMesh::EvaluatePS(NekDouble s)
int a = 0;
int b = 0;
ASSERTL1(!(s < 0) && !(s > m_curvelength),"s out of bounds");
ASSERTL1(!(s < 0) && !(s > m_curvelength), "s out of bounds");
if (s == 0)
{
......@@ -309,44 +342,30 @@ void CurveMesh::GetSampleFunction()
Array<OneD, NekDouble> loc = m_cadcurve->P(t);
/*NekDouble ts =
m_bl.Evaluate(m_blID, loc[0], loc[1], loc[2], 0.0);
bool found = false;
if (ts > 0.0)
// if inside the BL, dsti[0] set to the BL thickness, i.e. the offset
if (m_endoffset.count(0))
{
Array<OneD, NekDouble> N = m_cadcurve->N(t);
Array<OneD, NekDouble> Nwrt = m_cadcurve->NormalWRT(t, 0);
if(N[0]*N[0] + N[1]*N[1] + N[2]*N[2] < 1e-6)
if (dsti[1] < m_endoffset[0])
{
dsti[0] = m_mesh->m_octree->Query(loc);
dsti[0] = m_endoffset[0];
found = true;
}
else if ( N[0]*Nwrt[0] + N[1]*Nwrt[1] + N[2]*Nwrt[2] > 0)
{
//concave
dsti[0] = m_mesh->m_octree->Query(loc);
}
else
}
if (m_endoffset.count(1) && !found)
{
if (dsti[1] > m_curvelength - m_endoffset[1])
{
NekDouble R = 1.0 / m_cadcurve->Curvature(t);
if(R > 2.0*t)
{
R = 2.0*t;
}
Array<OneD, NekDouble> tloc(3);
tloc[0] = loc[0] + ts * Nwrt[0];
tloc[1] = loc[1] + ts * Nwrt[1];
tloc[2] = loc[2] + ts * Nwrt[2];
NekDouble d = m_mesh->m_octree->Query(tloc);
dsti[0] = d * R / (R + ts);
dsti[0] = m_endoffset[1];
found = true;
}
}
else
{*/
// else, dsti[0] is found from the octree
if (!found)
{
dsti[0] = m_mesh->m_octree->Query(loc);
//}
}
dsti[2] = t;
......@@ -356,7 +375,7 @@ void CurveMesh::GetSampleFunction()
void CurveMesh::PeriodicOverwrite(CurveMeshSharedPtr from)
{
//clear current mesh points and remove edges from edgeset
// clear current mesh points and remove edges from edgeset
m_meshpoints.clear();
for (int i = 0; i < m_meshedges.size(); i++)
{
......@@ -382,8 +401,8 @@ void CurveMesh::PeriodicOverwrite(CurveMeshSharedPtr from)
for (int i = 1; i < nodes.size() - 1; i++)
{
Array<OneD, NekDouble> loc = nodes[i]->GetLoc();
NodeSharedPtr nn = NodeSharedPtr(new Node(
m_mesh->m_numNodes++, loc[0] + T[0], loc[1] + T[1], 0.0));
NodeSharedPtr nn = NodeSharedPtr(
new Node(m_mesh->m_numNodes++, loc[0] + T[0], loc[1] + T[1], 0.0));
for (int j = 0; j < surfs.size(); j++)
{
......@@ -406,7 +425,7 @@ void CurveMesh::PeriodicOverwrite(CurveMeshSharedPtr from)
m_meshpoints.insert(m_meshpoints.begin(), verts[0]->GetNode());
m_meshpoints.push_back(verts[1]->GetNode());
//dont need to realign cad for vertices
// dont need to realign cad for vertices
// make edges and add them to the edgeset for the face mesher to use
for (int i = 0; i < m_meshpoints.size() - 1; i++)
......@@ -418,6 +437,5 @@ void CurveMesh::PeriodicOverwrite(CurveMeshSharedPtr from)
m_meshedges.push_back(e);
}
}
}
}
......@@ -136,6 +136,11 @@ public:
return m_id;
}
void SetOffset(unsigned i, NekDouble offset)
{
m_endoffset[i] = offset;
}
private:
/**
* @brief get node spacing sampling function
......@@ -188,6 +193,8 @@ private:
std::vector<NodeSharedPtr> m_meshpoints;
LibUtilities::AnalyticExpressionEvaluator m_bl;
int m_blID;
/// offset of second point at each end
std::map<unsigned, NekDouble> m_endoffset;
};
}
......
......@@ -3,7 +3,7 @@
<INFORMATION>
<I PROPERTY="CADFile" VALUE="2d-cad.STEP" />
<I PROPERTY="MeshType" VALUE="2D" />
<I PROPERTY="MeshType" VALUE="2DBndLayer" />
</INFORMATION>
<PARAMETERS>
......@@ -13,8 +13,8 @@
<P PARAM="Order" VALUE="4" />
<P PARAM="BndLayerSurfaces" VALUE="5-6" />
<P PARAM="BndLayerThickness" VALUE="0.02" />
<P PARAM="BndLayerSurfaces" VALUE="2,4" />
<P PARAM="BndLayerThickness" VALUE="0.1" />
<P PARAM="BndLayerLayers" VALUE="6" />
<P PARAM="BndLayerProgression" VALUE="1.5" />
</PARAMETERS>
......
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