Commit 5868d5be authored by Michael Turner's avatar Michael Turner
Browse files

added curve edge splitting high-order awareness and fixed issue with cad information propergation

parent 8873d0a4
......@@ -256,23 +256,35 @@ void CADSurf::Test(Array<OneD, NekDouble> uv)
if(uv[0] < m_occSurface.FirstUParameter())
{
error << "U(" << uv[0] << ") is less than Umin(" << m_occSurface.FirstUParameter() << ")";
passed = false;
if(fabs(uv[0]-m_occSurface.FirstUParameter()) > 1E-8)
{
error << "U(" << uv[0] << ") is less than Umin(" << m_occSurface.FirstUParameter() << ")";
passed = false;
}
}
else if(uv[0] > m_occSurface.LastUParameter())
{
error << "U(" << uv[0] << ") is greater than Umax(" << m_occSurface.LastUParameter() << ")";
passed = false;
if(fabs(uv[0]-m_occSurface.LastUParameter()) > 1E-8)
{
error << "U(" << uv[0] << ") is greater than Umax(" << m_occSurface.LastUParameter() << ")";
passed = false;
}
}
else if(uv[1] < m_occSurface.FirstVParameter())
{
error << "V(" << uv[1] << ") is less than Vmin(" << m_occSurface.FirstVParameter() << ")";
passed = false;
if(fabs(uv[1]-m_occSurface.FirstVParameter()) > 1E-8)
{
error << "V(" << uv[1] << ") is less than Vmin(" << m_occSurface.FirstVParameter() << ")";
passed = false;
}
}
else if(uv[1] > m_occSurface.LastVParameter())
{
error << "V(" << uv[1] << ") is greater than Vmax(" << m_occSurface.LastVParameter() << ")";
passed = false;
if(fabs(uv[1]-m_occSurface.LastVParameter()) > 1E-8)
{
error << "V(" << uv[1] << ") is greater than Vmax(" << m_occSurface.LastVParameter() << ")";
passed = false;
}
}
ASSERTL0(passed, error.str());
......
......@@ -138,7 +138,7 @@ namespace MeshUtils
vector<pair<ElementSharedPtr, int> > m_elLink;
/// id of cad curve which edge lies on
int CADCurveID;
int CADSurfID;
std::vector<int> CADSurfID;
private:
SpatialDomains::SegGeomSharedPtr m_geom;
......
......@@ -164,6 +164,51 @@ void CurveMesh::Mesh()
}
}
//post process the curve mesh to analyse for bad segments based on high-order normals and split if needed
int ct = 1;
while(ct > 0)
{
ct = 0;
for(int i = 0; i < m_meshpoints.size() - 1; i++)
{
bool split = false;
for(int j = 0; j < 2; j++)
{
Array<OneD, NekDouble> uv1, uv2;
uv1 = m_meshpoints[i]->GetCADSurf(s[j]->GetId());
uv2 = m_meshpoints[i+1]->GetCADSurf(s[j]->GetId());
Array<OneD, NekDouble> N1, N2;
N1 = s[j]->N(uv1);
N2 = s[j]->N(uv2);
NekDouble dot = N1[0]*N2[0] + N1[1]*N2[1] + N1[2]*N2[2];
if(acos(dot) > 3.142/2.0-0.1)
{
split = true;
}
}
if(split)
{
cout << "split" << endl;
ct++;
NekDouble t1, t2;
t1 = m_meshpoints[i]->GetCADCurve(m_id);
t2 = m_meshpoints[i+1]->GetCADCurve(m_id);
NekDouble tn = (t1 + t2)/2.0;
Array<OneD, NekDouble> loc = m_cadcurve->P(tn);
NodeSharedPtr nn = boost::shared_ptr<Node>(new Node(0,loc[0],loc[1],loc[2]));
nn->SetCADCurve(m_id, tn);
for(int j = 0; j < 2; j++)
{
Array<OneD, NekDouble> uv = s[j]->locuv(loc);
nn->SetCADSurf(s[j]->GetId(), uv);
}
m_meshpoints.insert(m_meshpoints.begin() + i+1, nn);
break;
}
}
}
}
......
......@@ -131,6 +131,11 @@ void FaceMesh::Mesh()
ElementSharedPtr E = GetElementFactory().CreateInstance(
LibUtilities::eTriangle, conf, m_connec[i], tags);
E->SetCADSurf(m_id);
vector<EdgeSharedPtr> edgs = E->GetEdgeList();
for(int j = 0; j < edgs.size(); j++)
{
edgs[j]->CADSurfID.push_back(m_id);
}
m_mesh->m_element[m_mesh->m_expDim].push_back(E);
}
}
......
......@@ -574,6 +574,7 @@ void SurfaceMesh::Optimise()
}
}
newe->CADSurfID.push_back(tri1->GetCADSurf());
m_mesh->m_edgeSet.insert(newe);
m_mesh->m_element[2][id1] = ntri1;
......
......@@ -199,6 +199,7 @@ void SurfaceMesh::HOSurf()
}
FaceSharedPtr f = m_mesh->m_element[2][i]->GetFaceLink();
vector<EdgeSharedPtr> edgelist = m_mesh->m_element[2][i]->GetEdgeList();
int surf = m_mesh->m_element[2][i]->GetTagList()[0];
vector<EdgeSharedPtr> egs = f->m_edgeList;
......@@ -214,8 +215,24 @@ void SurfaceMesh::HOSurf()
//insertion was sucsesful create high-order info
EdgeSharedPtr e = *testIns.first;
//the edges in the element are different to those in the face
//the cad information is stored in the element edges which are not
//in the m_mesh->m_edgeSet group.
//need to link them together and copy the cad information to be
//able to identify how to make it high-order
for(int k = 0; k < edgelist.size(); k++)
{
if(edgelist[k] == e)
{
e->CADSurfID = edgelist[k]->CADSurfID;
}
}
ASSERTL0(e->CADSurfID.size() == 1 || e->CADSurfID.size() == 2,
"incorrect cad surfs");
//figure out surfaces and curves this edge is on
if(e->m_n1->GetListCADCurve().size() > 0 && e->m_n2->GetListCADCurve().size() > 0)
if(e->CADSurfID.size()==2) //shoul be on curve
{
vector<int> clist1, clist2, c;
clist1 = e->m_n1->GetListCADCurve();
......@@ -246,16 +263,6 @@ void SurfaceMesh::HOSurf()
e->CADCurveID = c[0];
}
}
else
{
e->CADCurveID = -1;
e->CADSurfID = m_mesh->m_element[2][i]->GetTagList()[0];
}
}
else
{
e->CADCurveID = -1;
e->CADSurfID = m_mesh->m_element[2][i]->GetTagList()[0];
}
//if we are here the edge needs to be high-ordered
......@@ -298,8 +305,9 @@ void SurfaceMesh::HOSurf()
}
else
{
ASSERTL0(e->CADSurfID.size() == 1, "should only be one");
//edge is on surface and needs 2d optimisation
int sid = e->CADSurfID;
int sid = e->CADSurfID[0];
CADSurfSharedPtr s = m_cad->GetSurf(sid);
Array<OneD, NekDouble> uvb,uve;
uvb = e->m_n1->GetCADSurf(sid);
......
......@@ -164,6 +164,18 @@ void InputCAD::Process()
m_surfacemesh->HOSurf();
for(int i = 0; i < m_mesh->m_element[2].size(); i++)
{
FaceSharedPtr f = m_mesh->m_element[2][i]->GetFaceLink();
SpatialDomains::GeometrySharedPtr geom = f->GetGeom(m_mesh->m_spaceDim);
SpatialDomains::GeomFactorsSharedPtr gfac = geom->GetGeomFactors();
if(!gfac->IsValid())
{
cout << "warning: invalid face curavture in boundary element" << endl;
}
}
/*
//look over all surface elements
......
......@@ -204,6 +204,29 @@ void Module::ProcessEdges(bool ReprocessEdges)
}
}
if(ed->CADSurfID.size() > 0)
{
//sort out cad information
list<int> surfs;
for(int k = 0; k < ed->CADSurfID.size(); k++)
{
surfs.push_back(ed->CADSurfID[k]);
}
for(int k = 0; k < e2->CADSurfID.size(); k++)
{
surfs.push_back(e2->CADSurfID[k]);
}
surfs.sort(); surfs.unique();
e2->CADSurfID.clear();
list<int>::iterator it;
for(it = surfs.begin(); it != surfs.end(); it++)
{
e2->CADSurfID.push_back(*it);
}
ASSERTL0(e2->CADSurfID.size() == 1 || e2->CADSurfID.size() == 2,
"incorrect cad surfs");
}
// Update edge to element map.
e2->m_elLink.push_back(
pair<ElementSharedPtr,int>(elmt[i],j));
......
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