Commit 4de7bbff authored by Julian Marcon's avatar Julian Marcon
Browse files

Added non-O BL mesh support in 2D.

parent 09fc8931
......@@ -85,6 +85,7 @@ v4.4.0
- Add new two-dimensional mesher from NACA code or step file (!720)
- Add periodic boundary condition meshing in 2D (!733)
- Automatic peralign call if periodic boundary conditions present (!733)
- Support non-O BL meshing in 2D (!733)
- Fix inverted boundary layer in 2D (!736)
- More sensible element sizing with boudnary layers in 2D (!736)
- Change variable names in mcf file to make more sense (!736)
......
......@@ -121,35 +121,85 @@ void Generator2D::Process()
}
}
if(m_config["blcurves"].beenSet)
if (m_config["blcurves"].beenSet)
{
ParseUtils::GenerateSeqVector(m_config["blcurves"].as<string>().c_str(),
m_blCurves);
m_thickness_ID =
m_thickness.DefineFunction("x y z", m_config["blthick"].as<string>());
m_thickness_ID = m_thickness.DefineFunction(
"x y z", m_config["blthick"].as<string>());
}
// linear mesh all curves
for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); i++)
// linear mesh all curves with a BL
for (int i = 0; i < m_blCurves.size(); i++)
{
ASSERTL0(m_blCurves[i] <= m_mesh->m_cad->GetNumCurve(),
"Invalid curve defined with boundary layer");
if (m_mesh->m_verbose)
{
LibUtilities::PrintProgressbar(i, m_mesh->m_cad->GetNumCurve(),
LibUtilities::PrintProgressbar(i + 1, m_mesh->m_cad->GetNumCurve(),
"Curve progress");
}
m_curvemeshes[m_blCurves[i]] =
MemoryManager<CurveMesh>::AllocateSharedPtr(
m_blCurves[i], m_mesh, m_config["blthick"].as<string>());
m_curvemeshes[m_blCurves[i]]->Mesh();
}
LibUtilities::AnalyticExpressionEvaluator bl;
int blID = bl.DefineFunction("x y z", m_config["blthick"].as<string>());
// check curves with adjacent BLs
for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); i++)
{
vector<unsigned int>::iterator f =
find(m_blCurves.begin(), m_blCurves.end(), i);
if (f == m_blCurves.end())
if (f != m_blCurves.end())
{
m_curvemeshes[i] =
MemoryManager<CurveMesh>::AllocateSharedPtr(i, m_mesh);
continue;
}
else
vector<CADVertSharedPtr> verts =
m_mesh->m_cad->GetCurve(i)->GetVertex();
vector<NekDouble> offset(2);
offset[0] = 0.0;
offset[1] = 0.0;
for (int j = 0; j < 2; ++j)
{
m_curvemeshes[i] = MemoryManager<CurveMesh>::AllocateSharedPtr(
i, m_mesh, m_config["blthick"].as<string>());
NodeSharedPtr node = verts[j]->GetNode();
if (node->GetNumCadCurve())
{
Array<OneD, NekDouble> loc = node->GetLoc();
offset[j] = bl.Evaluate(blID, loc[0], loc[1], loc[2], 0.0);
}
}
m_curvemeshes[i] =
MemoryManager<CurveMesh>::AllocateSharedPtr(i, m_mesh, offset);
}
// linear mesh all other curves
for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); i++)
{
if (m_mesh->m_verbose)
{
LibUtilities::PrintProgressbar(i + m_blCurves.size(),
m_mesh->m_cad->GetNumCurve(),
"Curve progress");
}
vector<unsigned int>::iterator f =
find(m_blCurves.begin(), m_blCurves.end(), i);
if (f != m_blCurves.end())
{
continue;
}
m_curvemeshes[i]->Mesh();
......@@ -446,9 +496,30 @@ void Generator2D::MakeBL(int faceid, vector<EdgeLoopSharedPtr> e)
map<NodeSharedPtr, vector<EdgeSharedPtr> >::iterator it;
for (it = m_nodesToEdge.begin(); it != m_nodesToEdge.end(); it++)
{
ASSERTL0(it->second.size() > 0 && it->second.size() < 3,
"weirdness, most likely bl_surfs are incorrect");
if (it->second.size() < 2)
{
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);
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];
......@@ -525,6 +596,34 @@ void Generator2D::MakeBL(int faceid, vector<EdgeLoopSharedPtr> e)
m_mesh->m_element[2].push_back(E);
}
}
for (map<int, CurveMeshSharedPtr>::iterator ic = m_curvemeshes.begin();
ic != m_curvemeshes.end(); ++ic)
{
vector<unsigned>::iterator ib =
find(m_blCurves.begin(), m_blCurves.end(), ic->first);
if (ib != m_blCurves.end())
{
continue;
}
vector<NekDouble> bloffset = ic->second->GetBLOffset();
vector<NodeSharedPtr> nodes = ic->second->GetMeshPoints();
if (bloffset[0] != 0.0)
{
nodes.erase(nodes.begin());
}
if (bloffset[1] != 0.0)
{
nodes.erase(nodes.end() - 1);
}
ic->second = MemoryManager<CurveMesh>::AllocateSharedPtr(ic->first,
m_mesh, nodes);
}
}
void Generator2D::Report()
......
......@@ -48,9 +48,15 @@ void CurveMesh::Mesh()
m_bounds = m_cadcurve->Bounds();
m_curvelength = m_cadcurve->GetTotLength();
m_numSamplePoints =
int(m_curvelength / m_mesh->m_octree->GetMinDelta()) + 5;
ds = m_curvelength / (m_numSamplePoints - 1);
ASSERTL0(m_curvelength > m_bloffset[0] + m_bloffset[1],
"Boundary layers too thick for adjacent curve");
m_numSamplePoints = int((m_curvelength - m_bloffset[0] - m_bloffset[1]) /
m_mesh->m_octree->GetMinDelta()) +
5;
ds = (m_curvelength - m_bloffset[0] - m_bloffset[1]) /
(m_numSamplePoints - 1);
GetSampleFunction();
......@@ -63,28 +69,24 @@ void CurveMesh::Mesh()
Ne = round(Ae);
if (Ne + 1 < 2)
meshsvalue.clear();
meshsvalue.push_back(0.0);
if (m_bloffset[0] > 0.0)
{
meshsvalue.resize(2);
meshsvalue[0] = 0.0;
meshsvalue[1] = m_curvelength;
Ne = 1;
meshsvalue.push_back(m_bloffset[0]);
}
else
{
if (Ne > 1)
{
GetPhiFunction();
meshsvalue.resize(Ne + 1);
meshsvalue[0] = 0.0;
meshsvalue[Ne] = m_curvelength;
for (int i = 1; i < Ne; i++)
{
int iterationcounter = 0;
bool iterate = true;
int k = i;
NekDouble ski = meshsvalue[i - 1];
NekDouble ski = meshsvalue.back();
NekDouble lastSki;
while (iterate)
{
......@@ -100,10 +102,19 @@ void CurveMesh::Mesh()
ASSERTL0(iterationcounter < 1000000, "iteration failed");
}
meshsvalue[i] = ski;
meshsvalue.push_back(ski);
}
}
if (m_bloffset[1] > 0.0)
{
meshsvalue.push_back(m_bloffset[1]);
}
meshsvalue.push_back(m_curvelength);
Ne = meshsvalue.size() - 1;
NekDouble t;
Array<OneD, NekDouble> loc;
......@@ -205,11 +216,11 @@ NekDouble CurveMesh::EvaluateDS(NekDouble s)
int a = 0;
int b = 0;
if (s == 0)
if (s <= m_bloffset[0])
{
return m_dst[0][0];
}
else if (s == m_curvelength)
else if (s >= m_curvelength - m_bloffset[1])
{
return m_dst[m_numSamplePoints - 1][0];
}
......@@ -285,26 +296,19 @@ NekDouble CurveMesh::EvaluatePS(NekDouble s)
void CurveMesh::GetSampleFunction()
{
m_dst.resize(m_numSamplePoints);
Array<OneD, NekDouble> loc = m_cadcurve->P(m_bounds[0]);
Array<OneD, NekDouble> loc(3);
vector<NekDouble> dsti(3);
vector<NekDouble> dsti;
dsti.resize(3);
dsti[0] = m_mesh->m_octree->Query(loc);
dsti[1] = 0.0;
dsti[2] = m_bounds[0];
m_dst[0] = dsti;
for (int i = 1; i < m_numSamplePoints; i++)
for (int i = 0; i < m_numSamplePoints; i++)
{
dsti[1] = i * ds;
NekDouble t = m_cadcurve->tAtArcLength(dsti[1]);
dsti[1] = i * ds + m_bloffset[0];
NekDouble t = (i > 0 || m_bloffset[0] > 0.0)
? m_cadcurve->tAtArcLength(dsti[1])
: m_bounds[0];
loc = m_cadcurve->P(t);
NekDouble ts =
m_bl.Evaluate(m_blID, loc[0], loc[1], loc[2], 0.0);
NekDouble ts = m_bl.Evaluate(m_blID, loc[0], loc[1], loc[2], 0.0);
if (ts > 0.0)
{
......
......@@ -66,6 +66,10 @@ public:
{
m_blID = m_bl.DefineFunction("x y z", "0.0");
m_cadcurve = m_mesh->m_cad->GetCurve(m_id);
m_bloffset.resize(2);
m_bloffset[0] = 0.0;
m_bloffset[1] = 0.0;
}
/**
......@@ -92,6 +96,20 @@ public:
{
m_blID = m_bl.DefineFunction("x y z", expr);
m_cadcurve = m_mesh->m_cad->GetCurve(m_id);
m_bloffset.resize(2);
m_bloffset[0] = 0.0;
m_bloffset[1] = 0.0;
}
/**
* @brief alternative constructor with adjacent boundary layer curves
*/
CurveMesh(int id, MeshSharedPtr m, std::vector<NekDouble> bloffset)
: m_id(id), m_mesh(m), m_bloffset(bloffset)
{
m_blID = m_bl.DefineFunction("x y z", "0.0");
m_cadcurve = m_mesh->m_cad->GetCurve(m_id);
}
/**
......@@ -152,6 +170,14 @@ public:
return m_cadcurve;
}
/**
* @brief get the boundary layer offsets
*/
std::vector<NekDouble> GetBLOffset()
{
return m_bloffset;
}
private:
/**
* @brief get node spacing sampling function
......@@ -209,6 +235,8 @@ private:
std::vector<NodeSharedPtr> m_meshpoints;
LibUtilities::AnalyticExpressionEvaluator m_bl;
int m_blID;
/// boundary layer offset on vertices
std::vector<NekDouble> m_bloffset;
};
typedef boost::shared_ptr<CurveMesh> CurveMeshSharedPtr;
......
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