Commit 7fd20069 authored by Michael Turner's avatar Michael Turner
Browse files

core concept for 2D walls

parent 418ab898
......@@ -183,8 +183,10 @@ INCLUDE (ThirdPartyQT4)
INCLUDE (ThirdPartySMV)
INCLUDE (ThirdPartyTriangle)
INCLUDE (ThirdPartyTetGen)
INCLUDE (FindANN)
INCLUDE (ThirdPartyCCM)
INCLUDE (Doxygen)
IF( NEKTAR_USE_MKL )
......
......@@ -82,6 +82,7 @@ TARGET_LINK_LIBRARIES(NekMeshUtils LINK_PUBLIC SolverUtils)
IF(NEKTAR_USE_MESH)
TARGET_LINK_LIBRARIES(NekMeshUtils LINK_PRIVATE ${TRIANGLE_LIBRARY})
TARGET_LINK_LIBRARIES(NekMeshUtils LINK_PRIVATE ${TETGEN_LIBRARY})
TARGET_LINK_LIBRARIES(NekMeshUtils LINK_PRIVATE ${LIB_ANN})
ENDIF()
IF(NEKTAR_USE_MESH)
......
......@@ -40,6 +40,8 @@
#include <LocalRegions/MatrixKey.h>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <ANN/ANN.h>
using namespace std;
namespace Nektar
{
......@@ -143,8 +145,11 @@ void FaceMesh::Mesh()
void FaceMesh::MakeBL()
{
NekDouble min = m_octree->GetMinDelta();
//create intial boundary layer with thickness m_minDelta/2
for(int i = 1; i < orderedLoops.size(); i++)//dont do this to first loop
{
vector<blpair> blps;
for(int j = 0; j < orderedLoops[i].size(); j++)
{
//for each of the nodes make a new node which exists off the surface
......@@ -169,7 +174,7 @@ void FaceMesh::MakeBL()
Array<OneD, NekDouble> loc = orderedLoops[i][j]->GetLoc();
Array<OneD, NekDouble> tp(3);
for(int k = 0; k < 3; k++) tp[k] = m_bl*AN[k] + loc[k];
for(int k = 0; k < 3; k++) tp[k] = min/2.0*AN[k] + loc[k];
//project tp onto to surface to get new point
Array<OneD, NekDouble> uv(2);
m_cadsurf->ProjectTo(tp,uv);
......@@ -177,11 +182,120 @@ void FaceMesh::MakeBL()
tp[0],tp[1],tp[2]));
nn->SetCADSurf(m_id,m_cadsurf,uv);
blpairs.push_back(pair<NodeSharedPtr,NodeSharedPtr>(orderedLoops[i][j],nn));
blpair blp;
blp.first = orderedLoops[i][j];
blp.newn = nn;
blp.pos = nn;
blp.N = AN;
blp.stop = false;
blps.push_back(blp);
//place the new node into ordered loops for the surface triangulation to work With
orderedLoops[i][j] = nn;
}
blpairs.push_back(blps);
}
int numst = 9;
NekDouble dbl = (m_bl - min / 2) / numst;
for(int st = 0; st < numst; st++)
{
NekDouble blt = min /2.0 + st * dbl;
//advance the boundary layer point to their new possible locations
//then make ANN trees out of the possible locations
//then query the possible locations against the trees
//if the new position is not possible set flag to false and dont advance the point
//TODO add proximity search to look at points in own loop
for(int i = 0; i < blpairs.size(); i++)
{
for(int j = 0; j < blpairs[i].size(); j++)
{
if(blpairs[i][j].stop)
continue;
Array<OneD, NekDouble> loc = blpairs[i][j].first->GetLoc();
Array<OneD, NekDouble> tp(3);
for(int k = 0; k < 3; k++)
{
tp[k] = blt*blpairs[i][j].N[k] + loc[k];
}
//project tp onto to surface to get new point
Array<OneD, NekDouble> uv(2);
m_cadsurf->ProjectTo(tp,uv);
blpairs[i][j].pos->Move(tp,m_id,uv);
}
}
//make ann trees
vector<ANNkd_tree*> trees;
for(int i = 0; i < blpairs.size(); i++)
{
ANNpointArray dataPts;
ANNkd_tree* kdTree;
dataPts = annAllocPts(blpairs[i].size(), 3);
for(int j = 0; j < blpairs[i].size(); j++)
{
if(blpairs[i][j].stop)
{
dataPts[j][0] = blpairs[i][j].newn->m_x;
dataPts[j][1] = blpairs[i][j].newn->m_y;
dataPts[j][2] = blpairs[i][j].newn->m_z;
}
else
{
dataPts[j][0] = blpairs[i][j].pos->m_x;
dataPts[j][1] = blpairs[i][j].pos->m_y;
dataPts[j][2] = blpairs[i][j].pos->m_z;
}
}
kdTree = new ANNkd_tree(dataPts, blpairs[i].size(), 3);
trees.push_back(kdTree);
}
for(int i = 0; i < blpairs.size(); i++)
{
for(int j = 0; j < blpairs[i].size(); j++)
{
if(blpairs[i][j].stop)
continue;
ANNpoint queryPt;
ANNidxArray nnIdx;
ANNdistArray dists;
queryPt = annAllocPt(3);
nnIdx = new ANNidx[1];
dists = new ANNdist[1];
queryPt[0] = blpairs[i][j].pos->m_x;
queryPt[1] = blpairs[i][j].pos->m_y;
queryPt[2] = blpairs[i][j].pos->m_z;
for(int k = 0; k < trees.size(); k++)
{
if(i == k)
continue;
trees[k]->annkSearch(queryPt, 1, nnIdx, dists);
if(sqrt(dists[0]) < m_octree->Query(blpairs[i][j].first->GetLoc()))
{
blpairs[i][j].stop = true;
break;
}
else
{
blpairs[i][j].newn = blpairs[i][j].pos;
}
}
}
}
}
}
void FaceMesh::OptimiseLocalMesh()
......@@ -712,39 +826,43 @@ void FaceMesh::BuildLocalMesh()
*/
int tricomp = m_mesh->m_numcomp++;
//first build quads is bl surface
//first build quads if bl surface
if(m_makebl)
{
int quadcomp = m_mesh->m_numcomp++;
for(int i = 0; i < blpairs.size() - 1; i++)//this wont work for more than one sym plane loop
{
ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false);
vector<NodeSharedPtr> ns;
ns.push_back(blpairs[i].first);
ns.push_back(blpairs[i].second);
ns.push_back(blpairs[i+1].second);
ns.push_back(blpairs[i+1].first);
vector<int> tags;
tags.push_back(quadcomp);
ElementSharedPtr E = GetElementFactory().CreateInstance(
LibUtilities::eQuadrilateral, conf, ns, tags);
E->CADSurfId = m_id;
m_localElements.push_back(E);
}
ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false);
for(int i = 0; i < blpairs.size(); i++)
{
ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false);
vector<NodeSharedPtr> ns;
ns.push_back(blpairs.back().first);
ns.push_back(blpairs.back().second);
ns.push_back(blpairs[0].second);
ns.push_back(blpairs[0].first);
vector<int> tags;
tags.push_back(quadcomp);
ElementSharedPtr E = GetElementFactory().CreateInstance(
LibUtilities::eQuadrilateral, conf, ns, tags);
E->CADSurfId = m_id;
m_localElements.push_back(E);
for(int j = 0; j < blpairs[i].size() - 1; j++)
{
vector<NodeSharedPtr> ns;
ns.push_back(blpairs[i][j].first);
ns.push_back(blpairs[i][j].newn);
ns.push_back(blpairs[i][j+1].newn);
ns.push_back(blpairs[i][j+1].first);
vector<int> tags;
tags.push_back(quadcomp);
ElementSharedPtr E = GetElementFactory().CreateInstance(
LibUtilities::eQuadrilateral, conf, ns, tags);
E->CADSurfId = m_id;
m_localElements.push_back(E);
}
{
vector<NodeSharedPtr> ns;
ns.push_back(blpairs[i].back().first);
ns.push_back(blpairs[i].back().newn);
ns.push_back(blpairs[i][0].newn);
ns.push_back(blpairs[i][0].first);
vector<int> tags;
tags.push_back(quadcomp);
ElementSharedPtr E = GetElementFactory().CreateInstance(
LibUtilities::eQuadrilateral, conf, ns, tags);
E->CADSurfId = m_id;
m_localElements.push_back(E);
}
}
for(int i = 0; i < m_localElements.size(); i++)
{
vector<NodeSharedPtr> nods = m_localElements[i]->GetVertexList();
......
......@@ -47,6 +47,15 @@ namespace Nektar
namespace NekMeshUtils
{
struct blpair
{
NodeSharedPtr first;
NodeSharedPtr newn;
NodeSharedPtr pos;
Array<OneD, NekDouble> N;
bool stop;
};
/**
* @brief class for surface meshes on individual surfaces (paramter plane meshes)
*/
......@@ -172,8 +181,8 @@ private:
NekDouble m_bl;
/// should build boundary layer
bool m_makebl;
/// list of node links between node on loop and its corresponding interior node in quads
std::vector<std::pair<NodeSharedPtr, NodeSharedPtr> > blpairs;
/// list of bl information which follows the orderedLoops counting
std::vector<std::vector<blpair> > blpairs;
};
typedef boost::shared_ptr<FaceMesh> FaceMeshSharedPtr;
......
......@@ -78,6 +78,7 @@ void SurfaceMesh::Mesh()
//linear mesh all surfaces
for(int i = 1; i <= m_cad->GetNumSurf(); i++)
{
if(i > 7) return;
if(m_mesh->m_verbose)
{
LibUtilities::PrintProgressbar(i,m_cad->GetNumSurf(),
......
......@@ -213,7 +213,7 @@ void InputCAD::Process()
{
m_mesh->m_numcomp = 1; //just tets
}
//m_mesh->m_nummode = 2;
m_mesh->m_nummode = 2;
//create surface mesh
m_mesh->m_expDim--; //just to make it easier to surface mesh for now
......@@ -228,7 +228,23 @@ void InputCAD::Process()
ProcessElements ();
ProcessComposites();
m_surfacemesh->Report();
vector<ElementSharedPtr> el = m_mesh->m_element[m_mesh->m_expDim];
m_mesh->m_element[m_mesh->m_expDim].clear();
for(int i = 0; i < el.size(); i++)
{
if(el[i]->GetTagList()[0] == 6 || el[i]->GetTagList()[0] == 7)
{
m_mesh->m_element[m_mesh->m_expDim].push_back(el[i]);
}
}
ProcessVertices ();
ProcessEdges ();
ProcessFaces ();
ProcessElements ();
ProcessComposites();
/*m_surfacemesh->Report();
m_mesh->m_expDim = 3;
m_mesh->m_fields.push_back("u");
......@@ -272,7 +288,7 @@ void InputCAD::Process()
{
cout << endl;
cout << m_mesh->m_element[3].size() << endl;
}
}*/
}
}
......
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