Commit 1505cf65 authored by Michael Turner's avatar Michael Turner
Browse files

added small feature testing and angle anaylsis

parent 5c7c08bf
......@@ -135,6 +135,24 @@ Array<OneD, NekDouble> CADSurf::locuv(Array<OneD, NekDouble> p)
return uvr;
}
NekDouble CADSurf::DistanceTo(Array<OneD, NekDouble> p)
{
gp_Pnt loc(p[0] * 1000.0, p[1] * 1000.0, p[2] * 1000.0);
// alternative locuv methods
ShapeAnalysis_Surface sas(m_s);
sas.SetDomain(m_occSurface.FirstUParameter(),
m_occSurface.LastUParameter(),
m_occSurface.FirstVParameter(),
m_occSurface.LastVParameter());
gp_Pnt2d p2 = sas.ValueOfUV(loc,1e-7);
gp_Pnt p3 = sas.Value(p2);
return p3.Distance(loc);
}
bool CADSurf::IsPlane()
{
if(m_occSurface.GetType() == GeomAbs_Plane)
......
......@@ -166,6 +166,8 @@ public:
int GetId(){return m_ID;}
NekDouble DistanceTo(Array<OneD, NekDouble> p);
private:
void Test(Array<OneD, NekDouble> uv);
......
......@@ -61,6 +61,45 @@ void CADSystem::Report()
cout << "\tCAD has: " << m_surfs.size() << " surfaces." << endl;
}
void CADSystem::SmallFeatureAnalysis(NekDouble min)
{
vector<pair<int,NekDouble> >lens;
map<int, CADCurveSharedPtr>::iterator it;
for(it = m_curves.begin(); it != m_curves.end(); it++)
{
ASSERTL0(it->second->GetTotLength() > 1e-5,"curve too small for meshing");
if(it->second->GetTotLength() < 5.0 * min)
{
lens.push_back(pair<int,NekDouble>(
it->second->GetID(),it->second->GetTotLength()));
}
}
if(lens.size()>0)
{
cout << endl;
cout << "Small feature testing:" << endl;
for(int i = 0; i < lens.size(); i++)
{
cout << "\tCurve: " << lens[i].first << " has length: " << lens[i].second << endl;
}
if(lens.size() == 1)
cout << "\tThis curve's length is less than or close to minDelta" << endl;
else
cout << "\tThese curve's length are less than or close to minDelta" << endl;
cout << "\tIf these featurs are not planar the mesh generator may struggle ";
cout << "to produce a high-order mesh due to high curvature" << endl;
cout << "\tWould you like the Octree system to adjust accordingly, this may ";
cout << "compramise mesh smoothness (Y/N)" << endl;
string ans;
cin >> ans;
if(boost::iequals(ans,"Y"))
smallfeatreAdjust = true;
else
smallfeatreAdjust = false;
}
}
Array<OneD, NekDouble> CADSystem::GetBoundingBox()
{
Array<OneD, NekDouble> bound(6);
......
......@@ -151,6 +151,8 @@ public:
*/
bool InsideShape(Array<OneD, NekDouble> loc);
void SmallFeatureAnalysis(NekDouble min);
private:
/// Private function to add curve to CADSystem::m_verts.
void AddVert(int i, TopoDS_Shape in);
......@@ -168,6 +170,8 @@ private:
std::map<int, CADVertSharedPtr> m_verts;
/// occ master object
TopoDS_Shape shape;
bool smallfeatreAdjust;
};
typedef boost::shared_ptr<CADSystem> CADSystemSharedPtr;
......
......@@ -94,6 +94,94 @@ void SurfaceMesh::Mesh()
}
}
void SurfaceMesh::Metric()
{
NekDouble Angles[18] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
for(int i = 0; i < m_mesh->m_element[2].size(); i++)
{
vector<NodeSharedPtr> list = m_mesh->m_element[2][i]->GetVertexList();
vector<NekDouble> ang;
ang.push_back(list[0]->Angle(list[1],list[2]));
ang.push_back(list[1]->Angle(list[0],list[2]));
ang.push_back(list[2]->Angle(list[1],list[0]));
for(int j = 0; j < ang.size(); j++)
{
if(ang[j] < 10.0/180.0*3.142)
Angles[0]++;
else if(ang[j] < 20.0/180.0*3.142)
Angles[1]++;
else if(ang[j] < 30.0/180.0*3.142)
Angles[2]++;
else if(ang[j] < 40.0/180.0*3.142)
Angles[3]++;
else if(ang[j] < 50.0/180.0*3.142)
Angles[4]++;
else if(ang[j] < 60.0/180.0*3.142)
Angles[5]++;
else if(ang[j] < 70.0/180.0*3.142)
Angles[6]++;
else if(ang[j] < 80.0/180.0*3.142)
Angles[7]++;
else if(ang[j] < 90.0/180.0*3.142)
Angles[8]++;
else if(ang[j] < 100.0/180.0*3.142)
Angles[9]++;
else if(ang[j] < 110.0/180.0*3.142)
Angles[10]++;
else if(ang[j] < 120.0/180.0*3.142)
Angles[11]++;
else if(ang[j] < 130.0/180.0*3.142)
Angles[12]++;
else if(ang[j] < 140.0/180.0*3.142)
Angles[13]++;
else if(ang[j] < 150.0/180.0*3.142)
Angles[14]++;
else if(ang[j] < 160.0/180.0*3.142)
Angles[15]++;
else if(ang[j] < 170.0/180.0*3.142)
Angles[16]++;
else
Angles[17]++;
}
}
cout.unsetf(ios_base::floatfield);
if(m_mesh->m_verbose)
{
cout << endl << endl;
cout << "\tSurface Mesh Angle Histogram" << endl;
for(int i = 0; i < 18; i++)
{
cout << "\t\t" << i*10.0 << "\t<= ang < " << (i+1)*10.0 << ":\t" << Angles[i] << "\t";
for(int j = 0; j < ceil(Angles[i]/m_mesh->m_element[2].size()/3*80); j++)
cout << "=";
cout << endl;
}
}
/*//analyse t/R
EdgeSet::iterator it;
for(it = m_mesh->m_edgeSet.begin(); it != m_mesh->m_edgeSet.end(); it++)
{
EdgeSharedPtr e = *it;
int surf = e->CADSurf;
Array<OneD, NekDouble> loc1, loc2, locc(3);
loc1 = e->m_n1->GetCADSurf(surf);
loc2 = e->m_n2->GetCADSurf(surf);
locc[0] = (loc1[0] + loc2[0])/2.0;
locc[1] = (loc1[1] + loc2[1])/2.0;
locc[2] = (loc1[2] + loc2[2])/2.0;
CADSurfSharedPtr s = m_cad->GetSurf(surf);
NekDouble t = s->DistanceTo(locc);
}*/
}
//this mesh is valided that each egde is listed twice in the triangles
void SurfaceMesh::Validate()
{
......
......@@ -93,6 +93,8 @@ class SurfaceMesh
void HOAwareness();
void Metric();
private:
/**
......
......@@ -104,6 +104,8 @@ void InputCAD::Process()
m_cad->Report();
}
m_cad->SmallFeatureAnalysis(m_minDelta);
//create octree
OctreeSharedPtr m_octree = MemoryManager<Octree>::AllocateSharedPtr(m_cad,
m_mesh->m_verbose, m_minDelta,
......@@ -135,7 +137,11 @@ void InputCAD::Process()
m_surfacemesh->Optimise();
m_surfacemesh->HOAwareness();
//m_surfacemesh->HOAwareness();
m_surfacemesh->Metric();
exit(-1);
ClearElementLinks(); //mesh needs reprocessing to clean element and edge lists, easiest way to do it
ProcessVertices ();
......
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