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

tweak to tolerances

parent 7aeec892
......@@ -139,8 +139,10 @@ Array<OneD, NekDouble> BLMesh::GetNormal(vector<ElementSharedPtr> tris)
Np[2] /= mag;
NekDouble dot = 0.0;
int ct = 0;
while(fabs(dot - 1) > 1e-4)
{
ct++;
vector<NekDouble> a(N.size());
NekDouble aSum = 0.0;
for(int i = 0; i < N.size(); i++)
......@@ -180,6 +182,8 @@ Array<OneD, NekDouble> BLMesh::GetNormal(vector<ElementSharedPtr> tris)
Np[2] = 0.5* NpN[2] + (1.0-0.5)*Np[2];
dot = Np[0] * NpN[0] + Np[1] * NpN[1] + Np[2] * NpN[2];
ASSERTL0(ct < 10000,"failed to find normal");
}
return Np;
......@@ -211,6 +215,7 @@ void BLMesh::Mesh()
NodeSet::iterator it;
int ct = 0;
int failed = 0;
for(it = m_mesh->m_vertexSet.begin(); it != m_mesh->m_vertexSet.end(); it++, ct++)
{
//if (m_mesh->m_verbose)
......@@ -290,9 +295,6 @@ void BLMesh::Mesh()
if(Visability(g->second,bln.N) < 0.342)
{
cout << (*it)->m_id << "\t" << Visability(g->second,bln.N) << "\t";
cout << Visability(g->second,GetNormal(g->second)) << "\t";
Array<OneD, NekDouble> bestN = GetNormal(g->second);
NekDouble val = -1.0*numeric_limits<double>::max();
......@@ -330,15 +332,12 @@ void BLMesh::Mesh()
}
bln.N = bestN;
cout << Visability(g->second,bln.N);
if(Visability(g->second,bln.N) > 0.342)
{
cout << "\timproved" << endl;
}
else
if(Visability(g->second,bln.N) < 0.0)
{
cout << endl;
cout << "failed " << (*it)->m_x << " " << (*it)->m_y << " "
<< (*it)->m_z << endl;
failed++;
}
}
......@@ -367,6 +366,12 @@ void BLMesh::Mesh()
}
}
if(failed > 0)
{
cout << "some normals are not visible" << endl;
abort();
}
if (m_mesh->m_verbose)
{
cout << endl;
......@@ -459,70 +464,6 @@ void BLMesh::Mesh()
cout << endl;
}
//loop over all prisms, if invalid shrink until it is
//being careful to act on nodes which have already been shrunk
for(int i = 0; i < m_mesh->m_element[3].size(); i++)
{
if (m_mesh->m_verbose)
{
LibUtilities::PrintProgressbar(
i, m_mesh->m_element[3].size(), "Invalidity sweeping");
}
ElementSharedPtr el = m_mesh->m_element[3][i];
SpatialDomains::GeometrySharedPtr geom =
el->GetGeom(m_mesh->m_spaceDim);
SpatialDomains::GeomFactorsSharedPtr gfac =
geom->GetGeomFactors();
map<ElementSharedPtr, ElementSharedPtr>::iterator j = priToTri.find(el);
ASSERTL0(j != priToTri.end(), "not found");
vector<NodeSharedPtr> ns = j->second->GetVertexList();
while(!gfac->IsValid())
{
NekDouble maxbl = max(blData[ns[0]].bl, blData[ns[1]].bl);
maxbl = max(maxbl, blData[ns[2]].bl);
if(maxbl < 1E-6)
{
cout << "shrunk element too far, invalid mesh" << endl;
cout << ns[0]->m_id << endl;
cout << ns[1]->m_id << endl;
cout << ns[2]->m_id << endl;
break;
}
int ct = 0;
for(int j = 0; j < 3; j++)
{
map<NodeSharedPtr, blInfo>::iterator bli = blData.find(ns[j]);
ASSERTL0(bli != blData.end(), "not found");
if(bli->second.bl <= maxbl * 0.75)
{
ct++;
continue;
}
ASSERTL0(ct < 3,"skipped all 3");
bli->second.bl *= 0.75;
bli->second.pNode->m_x = ns[j]->m_x + bli->second.N[0] * bli->second.bl;
bli->second.pNode->m_y = ns[j]->m_y + bli->second.N[1] * bli->second.bl;
bli->second.pNode->m_z = ns[j]->m_z + bli->second.N[2] * bli->second.bl;
}
geom = el->GetGeom(m_mesh->m_spaceDim);
gfac = geom->GetGeomFactors();
}
}
if (m_mesh->m_verbose)
{
cout << endl;
}
//this is where it should do some clever collision dectecting and reduce the bl parameter
//all nodes in the vertex set are unique ordered and in surface elements
......@@ -675,6 +616,70 @@ void BLMesh::Mesh()
cout << endl;
}
//loop over all prisms, if invalid shrink until it is
//being careful to act on nodes which have already been shrunk
for(int i = 0; i < m_mesh->m_element[3].size(); i++)
{
if (m_mesh->m_verbose)
{
LibUtilities::PrintProgressbar(
i, m_mesh->m_element[3].size(), "Invalidity sweeping");
}
ElementSharedPtr el = m_mesh->m_element[3][i];
SpatialDomains::GeometrySharedPtr geom =
el->GetGeom(m_mesh->m_spaceDim);
SpatialDomains::GeomFactorsSharedPtr gfac =
geom->GetGeomFactors();
map<ElementSharedPtr, ElementSharedPtr>::iterator j = priToTri.find(el);
ASSERTL0(j != priToTri.end(), "not found");
vector<NodeSharedPtr> ns = j->second->GetVertexList();
while(!gfac->IsValid())
{
NekDouble maxbl = max(blData[ns[0]].bl, blData[ns[1]].bl);
maxbl = max(maxbl, blData[ns[2]].bl);
if(maxbl < 1E-6)
{
cout << "shrunk element too far, invalid mesh" << endl;
cout << ns[0]->m_id << endl;
cout << ns[1]->m_id << endl;
cout << ns[2]->m_id << endl;
break;
}
int ct = 0;
for(int j = 0; j < 3; j++)
{
map<NodeSharedPtr, blInfo>::iterator bli = blData.find(ns[j]);
ASSERTL0(bli != blData.end(), "not found");
if(bli->second.bl <= maxbl * 0.75)
{
ct++;
continue;
}
ASSERTL0(ct < 3,"skipped all 3");
bli->second.bl *= 0.75;
bli->second.pNode->m_x = ns[j]->m_x + bli->second.N[0] * bli->second.bl;
bli->second.pNode->m_y = ns[j]->m_y + bli->second.N[1] * bli->second.bl;
bli->second.pNode->m_z = ns[j]->m_z + bli->second.N[2] * bli->second.bl;
}
geom = el->GetGeom(m_mesh->m_spaceDim);
gfac = geom->GetGeomFactors();
}
}
if (m_mesh->m_verbose)
{
cout << endl;
}
//smoothing
//need to build a list of nodes to neigbours
map<ElementSharedPtr, ElementSharedPtr>::iterator eit;
......@@ -714,18 +719,6 @@ void BLMesh::Mesh()
}
}
/* vector<ElementSharedPtr> els = m_mesh->m_element[3];
m_mesh->m_element[3].clear();
for(int i = 0; i < els.size(); i++)
{
ElementSharedPtr tri = priToTri[els[i]];
if(tri->CADSurfId == 1 || tri->CADSurfId == 10)
{
m_mesh->m_element[3].push_back(els[i]);
}
}
*/
m_symSurfs = vector<int>(symSurfs.begin(), symSurfs.end());
//compile a map of the nodes needed for systemetry surfs
for(int i = 0; i < m_symSurfs.size(); i++)
......
......@@ -112,23 +112,19 @@ Array<OneD, NekDouble> CADSurf::locuv(Array<OneD, NekDouble> p)
uvr[1] < m_occSurface.FirstVParameter() ||
uvr[1] > m_occSurface.LastVParameter())
{
if (uvr[0] < m_occSurface.FirstUParameter() &&
fabs(m_occSurface.FirstUParameter() - uvr[0]) < 1E-6)
if (uvr[0] < m_occSurface.FirstUParameter() )
{
uvr[0] = m_occSurface.FirstUParameter();
}
else if (uvr[0] > m_occSurface.LastUParameter() &&
fabs(m_occSurface.LastUParameter() - uvr[0]) < 1E-6)
else if (uvr[0] > m_occSurface.LastUParameter() )
{
uvr[0] = m_occSurface.LastUParameter();
}
else if (uvr[1] < m_occSurface.FirstVParameter() &&
fabs(m_occSurface.FirstVParameter() - uvr[1]) < 1E-6)
else if (uvr[1] < m_occSurface.FirstVParameter())
{
uvr[1] = m_occSurface.FirstVParameter();
}
else if (uvr[1] > m_occSurface.LastVParameter() &&
fabs(m_occSurface.LastVParameter() - uvr[1]) < 1E-6)
else if (uvr[1] > m_occSurface.LastVParameter() )
{
uvr[1] = m_occSurface.LastVParameter();
}
......@@ -334,7 +330,7 @@ Array<OneD, NekDouble> CADSurf::D2(Array<OneD, NekDouble> uv)
void CADSurf::Test(Array<OneD, NekDouble> uv)
{
stringstream error;
/*stringstream error;
error << "Point not within parameter plane: ";
......@@ -342,7 +338,7 @@ void CADSurf::Test(Array<OneD, NekDouble> uv)
if (uv[0] < m_occSurface.FirstUParameter())
{
if (fabs(uv[0] - m_occSurface.FirstUParameter()) > 1E-8)
if (fabs(uv[0] - m_occSurface.FirstUParameter()) > 1E-6)
{
error << "U(" << uv[0] << ") is less than Umin("
<< m_occSurface.FirstUParameter() << ")";
......@@ -351,7 +347,7 @@ void CADSurf::Test(Array<OneD, NekDouble> uv)
}
else if (uv[0] > m_occSurface.LastUParameter())
{
if (fabs(uv[0] - m_occSurface.LastUParameter()) > 1E-8)
if (fabs(uv[0] - m_occSurface.LastUParameter()) > 1E-6)
{
error << "U(" << uv[0] << ") is greater than Umax("
<< m_occSurface.LastUParameter() << ")";
......@@ -360,7 +356,7 @@ void CADSurf::Test(Array<OneD, NekDouble> uv)
}
else if (uv[1] < m_occSurface.FirstVParameter())
{
if (fabs(uv[1] - m_occSurface.FirstVParameter()) > 1E-8)
if (fabs(uv[1] - m_occSurface.FirstVParameter()) > 1E-6)
{
error << "V(" << uv[1] << ") is less than Vmin("
<< m_occSurface.FirstVParameter() << ")";
......@@ -369,7 +365,7 @@ void CADSurf::Test(Array<OneD, NekDouble> uv)
}
else if (uv[1] > m_occSurface.LastVParameter())
{
if (fabs(uv[1] - m_occSurface.LastVParameter()) > 1E-8)
if (fabs(uv[1] - m_occSurface.LastVParameter()) > 1E-6)
{
error << "V(" << uv[1] << ") is greater than Vmax("
<< m_occSurface.LastVParameter() << ")";
......@@ -379,7 +375,7 @@ void CADSurf::Test(Array<OneD, NekDouble> uv)
error << " On Surface: " << GetId();
ASSERTL0(passed, error.str());
ASSERTL0(passed, error.str());*/
}
}
}
......@@ -73,7 +73,6 @@ Array<OneD, NekDouble> CADSystem::GetBoundingBox()
for (int i = 1; i <= m_curves.size(); i++)
{
gp_Pnt start, end;
CADCurveSharedPtr c = GetCurve(i);
Array<OneD, NekDouble> ends = c->GetMinMax();
......@@ -90,6 +89,62 @@ Array<OneD, NekDouble> CADSystem::GetBoundingBox()
return bound;
}
vector<int> CADSystem::GetBoundarySurfs()
{
vector<int> ret;
set<int> surfs;
vector<CADCurveSharedPtr> cs;
Array<OneD, NekDouble> bound = GetBoundingBox();
for (int i = 1; i <= m_curves.size(); i++)
{
CADCurveSharedPtr c = GetCurve(i);
Array<OneD, NekDouble> ends = c->GetMinMax();
if((fabs(bound[0] - ends[0]) < 1e-4 ||
fabs(bound[0] - ends[3]) < 1e-4 ||
fabs(bound[1] - ends[0]) < 1e-4 ||
fabs(bound[1] - ends[3]) < 1e-4)
&&
(fabs(bound[2] - ends[1]) < 1e-4 ||
fabs(bound[2] - ends[4]) < 1e-4 ||
fabs(bound[3] - ends[1]) < 1e-4 ||
fabs(bound[3] - ends[4]) < 1e-4)
&&
(fabs(bound[4] - ends[2]) < 1e-4 ||
fabs(bound[4] - ends[5]) < 1e-4 ||
fabs(bound[5] - ends[2]) < 1e-4 ||
fabs(bound[5] - ends[5]) < 1e-4) )
{
//curve touches on bounding box
cs.push_back(c);
}
}
for(int i = 0; i < cs.size(); i++)
{
vector<CADSurfSharedPtr> s = cs[i]->GetAdjSurf();
for(int j = 0; j < s.size(); j++)
{
surfs.insert(s[j]->GetId());
}
}
cout << endl;
set<int>::iterator it;
for(it = surfs.begin(); it != surfs.end(); it++)
{
ret.push_back(*it);
}
return ret;
}
bool CADSystem::LoadCAD()
{
if (!boost::filesystem::exists(m_name.c_str()))
......
......@@ -163,6 +163,8 @@ public:
*/
bool InsideShape(Array<OneD, NekDouble> loc);
std::vector<int> GetBoundarySurfs();
private:
/// Function to add curve to CADSystem::m_verts.
void AddVert(int i, TopoDS_Shape in);
......
......@@ -64,6 +64,11 @@ void TetGenInterface::InitialMesh(map<int, NodeSharedPtr> tgidton,
{
Array<OneD, NekDouble> loc = it->second->GetLoc();
if(it->first == 94127)
{
cout << loc[0] << " " << loc[1] << " " << loc[2] << endl;
}
surface.pointlist[it->first * 3 + 0] = loc[0];
surface.pointlist[it->first * 3 + 1] = loc[1];
surface.pointlist[it->first * 3 + 2] = loc[2];
......@@ -89,7 +94,7 @@ void TetGenInterface::InitialMesh(map<int, NodeSharedPtr> tgidton,
surface.facetmarkerlist[i] = 0;
}
tetrahedralize("pYzqQ", &surface, &output);
tetrahedralize("pYzQ", &surface, &output);
}
void TetGenInterface::GetNewPoints(int num,
......
......@@ -122,15 +122,6 @@ void TriangleInterface::Mesh(bool Quiet, bool Quality)
{
dt.triangulate("pYY");
}
// verify the mesh a bit
if (dt.out.numberofpoints - dt.out.numberofedges + dt.out.numberoftriangles !=
2 - m_centers.size())
{
cout << endl << "epc wrong" << endl;
cout << dt.out.numberofpoints - dt.out.numberofedges + dt.out.numberoftriangles
<< " " << m_centers.size() << " " << sid << endl;
}
}
void TriangleInterface::SetUp()
......
......@@ -197,7 +197,7 @@ void FaceMesh::QuadRemesh(map<NodeSharedPtr, NodeSharedPtr> nmap)
ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false);
vector<int> tags;
tags.push_back(m_id + 200);
tags.push_back(m_id + (over ? 2000 : 200));
ElementSharedPtr E = GetElementFactory().CreateInstance(
LibUtilities::eQuadrilateral, conf, ns, tags);
m_localElements.push_back(E);
......@@ -224,7 +224,7 @@ void FaceMesh::QuadRemesh(map<NodeSharedPtr, NodeSharedPtr> nmap)
ElmtConfig conf(LibUtilities::eQuadrilateral, 1, false, false);
vector<int> tags;
tags.push_back(m_id + 200);
tags.push_back(m_id + (over ? 2000 : 200));
ElementSharedPtr E = GetElementFactory().CreateInstance(
LibUtilities::eQuadrilateral, conf, ns, tags);;
m_localElements.push_back(E);
......@@ -429,11 +429,12 @@ void FaceMesh::DiagonalSwap()
NodeSet::iterator nit;
for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
{
if ((*nit)->GetNumCadCurve() == 0)
{
//this routine is broken and needs looking at
//if ((*nit)->GetNumCadCurve() == 0)
//{
// node is interior
idealConnec[(*nit)->m_id] = 6;
}
/*}
else
{
// need to identify the two other nodes on the boundary to find
......@@ -456,7 +457,7 @@ void FaceMesh::DiagonalSwap()
idealConnec[(*nit)->m_id] =
ceil((*nit)->Angle(ns[0], ns[1]) / 3.142 * 3) + 1;
}
}*/
}
for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
{
......@@ -806,7 +807,7 @@ void FaceMesh::BuildLocalMesh()
ElmtConfig conf(LibUtilities::eTriangle, 1, false, false);
vector<int> tags;
tags.push_back(m_id + 100);
tags.push_back(m_id + (over ? 1000 : 100));
ElementSharedPtr E = GetElementFactory().CreateInstance(
LibUtilities::eTriangle, conf, m_connec[i], tags);
E->CADSurfId = m_id;
......@@ -1089,8 +1090,8 @@ void FaceMesh::OrientateCurves()
N[1] = (n2info[0] - n1info[0]) / mag;
Array<OneD, NekDouble> P(2);
P[0] = (n1info[0] + n2info[0]) / 2.0 + 1e-6 * N[0];
P[1] = (n1info[1] + n2info[1]) / 2.0 + 1e-6 * N[1];
P[0] = (n1info[0] + n2info[0]) / 2.0 + 1e-8 * N[0];
P[1] = (n1info[1] + n2info[1]) / 2.0 + 1e-8 * N[1];
// now test to see if p is inside or outside the shape
// vector to the right
......
......@@ -62,9 +62,10 @@ public:
MeshSharedPtr m,
CADSurfSharedPtr cad,
OctreeSharedPtr oct,
const std::map<int, CurveMeshSharedPtr> &cmeshes)
const std::map<int, CurveMeshSharedPtr> &cmeshes,
bool ov)
: m_mesh(m), m_cadsurf(cad), m_octree(oct), m_curvemeshes(cmeshes),
m_id(id)
m_id(id), over(ov)
{
m_edgeloops = m_cadsurf->GetEdges();
......@@ -152,6 +153,8 @@ private:
EdgeSet m_localEdges;
/// local list of elements
std::vector<ElementSharedPtr> m_localElements;
bool over;
};
typedef boost::shared_ptr<FaceMesh> FaceMeshSharedPtr;
......
......@@ -75,6 +75,9 @@ void SurfaceMesh::Mesh()
// linear mesh all surfaces
for (int i = 1; i <= m_cad->GetNumSurf(); i++)
{
//if(i == 1200 || i == 2206)
// continue;
if (m_mesh->m_verbose)
{
LibUtilities::PrintProgressbar(
......@@ -83,7 +86,7 @@ void SurfaceMesh::Mesh()
m_facemeshes[i] =
MemoryManager<FaceMesh>::AllocateSharedPtr(i,m_mesh,
m_cad->GetSurf(i), m_octree, m_curvemeshes);
m_cad->GetSurf(i), m_octree, m_curvemeshes, m_cad->GetNumSurf() > 100);
m_facemeshes[i]->Mesh();
}
......
......@@ -397,7 +397,6 @@ void SurfaceMesh::HOSurf()
ti[k] = tb * (1.0 - gll[k]) / 2.0 +
te * (1.0 + gll[k]) / 2.0;
}
exit(-1);
break;
}
itct++;
......
......@@ -114,18 +114,6 @@ void InputCAD::Process()
m_writeoctree = pSession->GetSolverInfo("WriteOctree") == "TRUE";
}
vector<unsigned int> symsurfs;
vector<unsigned int> blsurfs;
if (m_makeBL)
{
string sym, bl;
bl = pSession->GetSolverInfo("BLSurfs");
ParseUtils::GenerateSeqVector(bl.c_str(), blsurfs);
sort(blsurfs.begin(), blsurfs.end());
ASSERTL0(blsurfs.size() > 0,
"No surfaces selected to make boundary layer on");
}
CADSystemSharedPtr m_cad =
MemoryManager<CADSystem>::AllocateSharedPtr(m_CADName);
......@@ -136,6 +124,37 @@ void InputCAD::Process()
ASSERTL0(m_cad->LoadCAD(), "Failed to load CAD");
vector<int> bs = m_cad->GetBoundarySurfs();
vector<unsigned int> symsurfs;
vector<unsigned int> blsurfs, blsurfst;
if (m_makeBL)
{
string sym, bl;
bl = pSession->GetSolverInfo("BLSurfs");
ParseUtils::GenerateSeqVector(bl.c_str(), blsurfst);
sort(blsurfst.begin(), blsurfst.end());
ASSERTL0(blsurfst.size() > 0,
"No surfaces selected to make boundary layer on");
for(int i = 0; i < blsurfst.size(); i++)
{
bool add = true;
for(int j = 0; j < bs.size(); j++)
{
if(bs[j] == blsurfst[i])
{
add = false;
break;
}
}
if(add)
{
blsurfs.push_back(blsurfst[i]);
}
}
}