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

improved angle based point insertion

parent a964cfd5
......@@ -40,8 +40,8 @@
#include <iomanip>
#include <SpatialDomains/PointGeom.h>
#include <NekMeshUtils/CADSystem/CADSystem.h>
#include <SpatialDomains/PointGeom.h>
namespace Nektar
{
......@@ -148,7 +148,7 @@ public:
NEKMESHUTILS_EXPORT NodeSharedPtr copy()
{
return boost::shared_ptr<Node>(new Node(m_id,m_x,m_y,m_z));
return boost::shared_ptr<Node>(new Node(m_id, m_x, m_y, m_z));
}
NEKMESHUTILS_EXPORT NekDouble abs2() const
......@@ -163,8 +163,7 @@ public:
NEKMESHUTILS_EXPORT Node curl(const Node &pSrc) const
{
return Node(m_id,
m_y * pSrc.m_z - m_z * pSrc.m_y,
return Node(m_id, m_y * pSrc.m_z - m_z * pSrc.m_y,
m_z * pSrc.m_x - m_x * pSrc.m_z,
m_x * pSrc.m_y - m_y * pSrc.m_x);
}
......@@ -228,7 +227,7 @@ public:
return an;
}
// fucntions for cad information
// fucntions for cad information
void SetCADCurve(int i, CADCurveSharedPtr c, NekDouble t)
{
......@@ -292,8 +291,8 @@ public:
std::vector<std::pair<int, CADSurfSharedPtr> > GetCADSurfs()
{
std::vector<std::pair<int, CADSurfSharedPtr> > lst;
std::map<int, std::pair<CADSurfSharedPtr, Array<OneD, NekDouble> > >::
iterator s;
std::map<int, std::pair<CADSurfSharedPtr,
Array<OneD, NekDouble> > >::iterator s;
for (s = CADSurfList.begin(); s != CADSurfList.end(); s++)
{
lst.push_back(
......@@ -324,44 +323,80 @@ public:
void MoveCurve(Array<OneD, NekDouble> l, int s, NekDouble t)
{
m_x = l[0];
m_y = l[1];
m_z = l[2];
m_x = l[0];
m_y = l[1];
m_z = l[2];
CADCurveSharedPtr cu = CADCurveList[s].first;
CADCurveList[s] =
std::pair<CADCurveSharedPtr, NekDouble >(cu, t);
CADCurveList[s] = std::pair<CADCurveSharedPtr, NekDouble>(cu, t);
}
std::vector<std::pair<int,std::string> > GetCADCurveInfoVector()
std::vector<std::pair<int, std::string> > GetCADCurveInfoVector()
{
std::vector<std::pair<int,std::string> > ret;
std::vector<std::pair<int, std::string> > ret;
std::map<int, std::pair<CADCurveSharedPtr, NekDouble> >::iterator c;
for (c = CADCurveList.begin(); c != CADCurveList.end(); c++)
{
std::stringstream ss;
ss << std::scientific << std::setprecision(8);
ss << c->second.second;
ret.push_back(std::pair<int,std::string>(c->first,ss.str()));
ret.push_back(std::pair<int, std::string>(c->first, ss.str()));
}
return ret;
}
std::vector<std::pair<int,std::string> > GetCADSurfInfoVector()
std::vector<std::pair<int, std::string> > GetCADSurfInfoVector()
{
std::vector<std::pair<int,std::string> > ret;
std::map<int, std::pair<CADSurfSharedPtr, Array<OneD, NekDouble> > >::
iterator s;
std::vector<std::pair<int, std::string> > ret;
std::map<int, std::pair<CADSurfSharedPtr,
Array<OneD, NekDouble> > >::iterator s;
for (s = CADSurfList.begin(); s != CADSurfList.end(); s++)
{
std::stringstream ss;
ss << std::scientific << std::setprecision(8);
ss << s->second.second[0] << " " << s->second.second[1];
ret.push_back(std::pair<int,std::string>(s->first,ss.str()));
ret.push_back(std::pair<int, std::string>(s->first, ss.str()));
}
return ret;
}
NekDouble Angle(Array<OneD, NekDouble> locA, Array<OneD, NekDouble> locB,
Array<OneD, NekDouble> N, int sid)
{
// calculates the angle between this node to a to this node to b
// Uses the CAD surface to orientate the angle
Array<OneD, NekDouble> A(3), B(3), CP(3);
A[0] = locA[0] - m_x;
A[1] = locA[1] - m_y;
A[2] = locA[2] - m_z;
B[0] = locB[0] - m_x;
B[1] = locB[1] - m_y;
B[2] = locB[2] - m_z;
CP[0] = A[1] * B[2] - A[2] * B[1];
CP[1] = -1.0 * (A[0] * B[2] - A[2] * B[0]);
CP[2] = A[0] * B[1] - A[1] * B[0];
NekDouble ang = sqrt(CP[0] * CP[0] + CP[1] * CP[1] + CP[2] * CP[2]);
ang /= sqrt(A[0] * A[0] + A[1] * A[1] + A[2] * A[2]);
ang /= sqrt(B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
std::map<int, std::pair<CADSurfSharedPtr,
Array<OneD, NekDouble> > >::iterator it;
it = CADSurfList.find(sid);
NekDouble dot = N[0] * CP[0] + N[1] * CP[1] + N[2] * CP[2];
ang = asin(ang);
if (dot < 0.0)
{
ang = 2.0 * M_PI - ang;
}
return ang;
}
/// ID of node.
int m_id;
/// X-coordinate.
......@@ -408,7 +443,6 @@ struct NodeHash : std::unary_function<NodeSharedPtr, std::size_t>
}
};
typedef boost::unordered_set<NodeSharedPtr, NodeHash> NodeSet;
}
}
......
......@@ -183,14 +183,14 @@ void FaceMesh::Mesh()
pplanemesh->Mesh();
pplanemesh->Extract(m_connec);
meshcounter++;
break;
// break;
}
// build a local version of the mesh (one set of triangles). this is done
// so edge connectivity infomration can be used for optimisation
BuildLocalMesh();
//OptimiseLocalMesh();
OptimiseLocalMesh();
// make new elements and add to list from list of nodes and connectivity
// from triangle removing unnesercary infomration from the elements
......@@ -904,22 +904,31 @@ bool FaceMesh::Validate()
int pointBefore = m_stienerpoints.size();
for (int i = 0; i < m_connec.size(); i++)
{
Array<OneD, NekDouble> r(3);
Array<OneD, NekDouble> r(3), a(3);
vector<Array<OneD, NekDouble> > info;
for(int j = 0; j < 3; j++)
{
info.push_back(m_connec[i][j]->GetCADSurfInfo(m_id));
}
r[0] = m_connec[i][0]->Distance(m_connec[i][1]);
r[1] = m_connec[i][1]->Distance(m_connec[i][2]);
r[2] = m_connec[i][2]->Distance(m_connec[i][0]);
a[0] = m_connec[i][0]->Angle(m_connec[i][1]->GetLoc(),
m_connec[i][2]->GetLoc(), m_cadsurf->N(info[0]), m_id);
a[1] = m_connec[i][1]->Angle(m_connec[i][2]->GetLoc(),
m_connec[i][0]->GetLoc(), m_cadsurf->N(info[1]), m_id);
a[2] = m_connec[i][2]->Angle(m_connec[i][0]->GetLoc(),
m_connec[i][1]->GetLoc(), m_cadsurf->N(info[2]), m_id);
NekDouble d1 = m_mesh->m_octree->Query(m_connec[i][0]->GetLoc());
NekDouble d2 = m_mesh->m_octree->Query(m_connec[i][1]->GetLoc());
NekDouble d3 = m_mesh->m_octree->Query(m_connec[i][2]->GetLoc());
vector<Array<OneD, NekDouble> > info;
for(int j = 0; j < 3; j++)
{
info.push_back(m_connec[i][j]->GetCADSurfInfo(m_id));
}
Array<OneD, NekDouble> uvc(2);
uvc[0] = (info[0][0] + info[1][0] + info[2][0]) / 3.0;
......@@ -931,29 +940,39 @@ bool FaceMesh::Validate()
NekDouble d = (d1 + d2 + d3 + d4) / 4.0;
vector<bool> valid(3);
valid[0] = r[0] < d * 1.41;
valid[1] = r[1] < d * 1.41;
valid[2] = r[2] < d * 1.41;
valid[0] = r[0] < d * 1.5;
valid[1] = r[1] < d * 1.5;
valid[2] = r[2] < d * 1.5;
vector<bool> angValid(3);
angValid[0] = a[0]/M_PI*180.0 > 20.0 && a[0]/M_PI*180.0 < 120.0;
angValid[1] = a[1]/M_PI*180.0 > 20.0 && a[1]/M_PI*180.0 < 120.0;
angValid[2] = a[2]/M_PI*180.0 > 20.0 && a[2]/M_PI*180.0 < 120.0;
int numValid = 0;
int numAngValid = 0;
for(int j = 0; j < 3; j++)
{
if(valid[j])
{
numValid++;
}
if(angValid[j])
{
numAngValid++;
}
}
cout << numValid << endl;
//if numvalid is zero no work to be done
/*if (numValid != 3)
{
AddNewPoint(uvc);
}*/
if(numValid != 3)
if(numValid != 3 || numAngValid != 3)
{
//break the bad edge
int a=0, b=0;
/*int a=0, b=0;
if(!valid[0])
{
a = 0;
......@@ -972,7 +991,7 @@ bool FaceMesh::Validate()
Array<OneD, NekDouble> uvn(2);
uvn[0] = (info[a][0] + info[b][0]) / 2.0;
uvn[1] = (info[a][1] + info[b][1]) / 2.0;
uvn[1] = (info[a][1] + info[b][1]) / 2.0;*/
AddNewPoint(uvc);
}
}
......
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