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

fix 2d wing case

parent 25a08b67
......@@ -113,7 +113,7 @@ void Generator2D::Process()
for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
{
//MakeBL(i);
MakeBL(i);
}
}
......@@ -217,8 +217,13 @@ void Generator2D::MakeBL(int faceid)
swap(p1, p2);
}
Array<OneD, NekDouble> n(2);
n[0] = p1[1] - p2[1];
n[1] = p2[0] - p1[0];
n[0] = p2[1] - p1[1];
n[1] = p1[0] - p2[0];
if(m_mesh->m_cad->GetSurf(faceid)->IsReversedNormal())
{
n[0] *= -1.0;
n[1] *= -1.0;
}
NekDouble mag = sqrt(n[0] * n[0] + n[1] * n[1]);
n[0] /= mag;
n[1] /= mag;
......
......@@ -56,7 +56,7 @@ void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
{
// this piece of code orientates the surface,
// it used to be face mesh but its easier to have it here
int np = 40;
int np = 20;
vector<vector<Array<OneD, NekDouble> > > loopt;
for (int i = 0; i < ein.size(); i++)
{
......@@ -93,12 +93,14 @@ void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
{
NekDouble area = 0.0;
bg::model::polygon<point_xy> polygon;
bg::model::polygon<point_xy, false, true> polygon;
vector<point_xy> points;
for(int j = 0; j < loopt[i].size(); j++)
{
points.push_back(point_xy(loopt[i][j][0], loopt[i][j][1]));
}
//boost requires for closed polygons (last point == first point)
points.push_back(point_xy(loopt[i][0][0], loopt[i][0][1]));
bg::assign_points(polygon,points);
......@@ -168,34 +170,33 @@ void CADSurf::OrientateEdges(CADSurfSharedPtr surf,
}
}
for(int i = 0; i < ein.size(); i++)
//only need center points for inner loops
for(int i = 1; i < ein.size(); i++)
{
Array<OneD, NekDouble> n1 = loopt[i][0];
Array<OneD, NekDouble> n2 = loopt[i][1];
if(ein[i]->edgeo[i] == CADSystem::eForwards)
{
swap(n1, n2);
}
Array<OneD, NekDouble> N(2);
NekDouble mag = sqrt((n1[0] - n2[0]) * (n1[0] - n2[0]) +
(n1[1] - n2[1]) * (n1[1] - n2[1]));
N[0] = -1.0 * (n2[1] - n1[1]) / mag;
N[1] = (n2[0] - n1[0]) / mag;
N[0] = (n2[1] - n1[1]) / mag;
N[1] = -1.0 * (n2[0] - n1[0]) / mag;
Array<OneD, NekDouble> P(2);
P[0] = (n1[0] + n2[0]) / 2.0 + 1e-6 * N[0];
P[1] = (n1[1] + n2[1]) / 2.0 + 1e-6 * N[1];
P[0] = (n1[0] + n2[0]) / 2.0 + N[0];
P[1] = (n1[1] + n2[1]) / 2.0 + N[1];
ein[i]->center = P;
bg::model::polygon<point_xy> polygon;
bg::model::polygon<point_xy, false, true> polygon;
vector<point_xy> points;
for(int j = 0; j < loopt[i].size(); j++)
{
points.push_back(point_xy(loopt[i][j][0], loopt[i][j][1]));
}
//boost requires for closed polygons (last point == first point)
points.push_back(point_xy(loopt[i][0][0], loopt[i][0][1]));
point_xy p(P[0],P[1]);
......
......@@ -60,7 +60,7 @@ public:
/**
* @brief Default constructor.
*/
CADSurf()
CADSurf() : m_correctNormal(true)
{
m_type = CADType::eSurf;
}
......@@ -157,7 +157,7 @@ public:
*/
bool IsReversedNormal()
{
return !m_correctNormal;
return m_correctNormal;
}
protected:
......
......@@ -48,7 +48,7 @@ std::string CADCurveOCE::key = GetCADCurveFactory().RegisterCreatorFunction(
NekDouble CADCurveOCE::tAtArcLength(NekDouble s)
{
NekDouble dt =
(m_occCurve.LastParameter() - m_occCurve.FirstParameter()) / (5000);
(m_occCurve.LastParameter() - m_occCurve.FirstParameter()) / (1000);
NekDouble t = m_occCurve.FirstParameter();
NekDouble len = 0.0;
......
......@@ -54,11 +54,6 @@ void CADSurfOCE::Initialise(int i, TopoDS_Shape in)
//reverse the face
if(in.Orientation() == 0)
{
BRepBuilderAPI_MakeFace nf;
m_s->VReverse();
nf.Init(m_s,true,1e-6);
in = nf.Face();
SetReverseNomral();
}
......@@ -92,8 +87,7 @@ Array<OneD, NekDouble> CADSurfOCE::locuv(Array<OneD, NekDouble> p)
GeomAPI_ProjectPointOnSurf projection(
loc, m_s, m_occSurface.FirstUParameter(), m_occSurface.LastUParameter(),
m_occSurface.FirstVParameter(), m_occSurface.LastVParameter(),
Extrema_ExtAlgo_Tree);
m_occSurface.FirstVParameter(), m_occSurface.LastVParameter());
Array<OneD, NekDouble> uvr(2);
if (projection.NbPoints() == 0)
......@@ -104,7 +98,7 @@ Array<OneD, NekDouble> CADSurfOCE::locuv(Array<OneD, NekDouble> p)
m_occSurface.FirstUParameter(), m_occSurface.LastUParameter(),
m_occSurface.FirstVParameter(), m_occSurface.LastVParameter());
gp_Pnt2d p2 = sas.ValueOfUV(loc, 1e-7);
gp_Pnt2d p2 = sas.ValueOfUV(loc, 1e-3);
uvr[0] = p2.X();
uvr[1] = p2.Y();
......
......@@ -165,7 +165,8 @@ bool CADSystemOCE::LoadCAD()
int e = mapOfEdges.FindIndex(edge);
edgeloop->edges.push_back(m_curves[e]);
edgeloop->edgeo.push_back(
exp.Orientation() == 0 ? eForwards : eBackwards);
exp.Orientation() == TopAbs_FORWARD ? eForwards
: eBackwards);
}
exp.Next();
......
......@@ -837,49 +837,52 @@ int Octree::CountElemt()
void Octree::CompileSourcePointList()
{
for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); i++)
if(m_mesh->m_cad->Is2D())
{
CADCurveSharedPtr curve = m_mesh->m_cad->GetCurve(i);
Array<OneD, NekDouble> bds = curve->GetBounds();
int samples = 100;
NekDouble dt = (bds[1] - bds[0]) / (samples + 1);
for (int j = 1; j < samples - 1; j++) // dont want first and last point
for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); i++)
{
NekDouble t = bds[0] + dt * j;
NekDouble C = curve->Curvature(t);
Array<OneD, NekDouble> loc = curve->P(t);
CADCurveSharedPtr curve = m_mesh->m_cad->GetCurve(i);
Array<OneD, NekDouble> bds = curve->GetBounds();
int samples = 100;
NekDouble dt = (bds[1] - bds[0]) / (samples + 1);
for (int j = 1; j < samples - 1; j++) // dont want first and last point
{
NekDouble t = bds[0] + dt * j;
NekDouble C = curve->Curvature(t);
vector<pair<CADSurfSharedPtr, CADSystem::Orientation> > ss =
curve->GetAdjSurf();
Array<OneD, NekDouble> uv = ss[0].first->locuv(loc);
Array<OneD, NekDouble> loc = curve->P(t);
if (C != 0.0)
{
NekDouble del = 2.0 * (1.0 / C) * sqrt(m_eps * (2.0 - m_eps));
vector<pair<CADSurfSharedPtr, CADSystem::Orientation> > ss =
curve->GetAdjSurf();
Array<OneD, NekDouble> uv = ss[0].first->locuv(loc);
if (del > m_maxDelta)
{
del = m_maxDelta;
}
if (del < m_minDelta)
if (C != 0.0)
{
del = m_minDelta;
}
NekDouble del = 2.0 * (1.0 / C) * sqrt(m_eps * (2.0 - m_eps));
CPointSharedPtr newCPoint =
MemoryManager<CPoint>::AllocateSharedPtr(
ss[0].first->GetId(), uv, loc, del);
if (del > m_maxDelta)
{
del = m_maxDelta;
}
if (del < m_minDelta)
{
del = m_minDelta;
}
m_SPList.push_back(newCPoint);
}
else
{
BPointSharedPtr newBPoint =
MemoryManager<BPoint>::AllocateSharedPtr(
ss[0].first->GetId(), uv, loc);
CPointSharedPtr newCPoint =
MemoryManager<CPoint>::AllocateSharedPtr(
ss[0].first->GetId(), uv, loc, del);
m_SPList.push_back(newCPoint);
}
else
{
BPointSharedPtr newBPoint =
MemoryManager<BPoint>::AllocateSharedPtr(
ss[0].first->GetId(), uv, loc);
m_SPList.push_back(newBPoint);
m_SPList.push_back(newBPoint);
}
}
}
}
......
......@@ -49,7 +49,7 @@ void CurveMesh::Mesh()
m_bounds = m_cadcurve->GetBounds();
m_curvelength = m_cadcurve->GetTotLength();
m_numSamplePoints =
int(m_curvelength / m_mesh->m_octree->GetMinDelta()) * 2;
int(m_curvelength / m_mesh->m_octree->GetMinDelta()) + 10;
ds = m_curvelength / (m_numSamplePoints - 1);
GetSampleFunction();
......
......@@ -373,7 +373,7 @@ void FaceMesh::Smoothing()
u0[1] += uj[1] / nodesystem.size();
}
Array<OneD, NekDouble> pu0 = m_cadsurf->P(u0);
/*Array<OneD, NekDouble> pu0 = m_cadsurf->P(u0);
NekDouble di = m_mesh->m_octree->Query(pu0);
Array<OneD, NekDouble> F(2, 0.0), dF(4, 0.0);
for (int i = 0; i < nodesystem.size(); i++)
......@@ -429,7 +429,7 @@ void FaceMesh::Smoothing()
dF[2] *= -1.0 / det;
u0[0] -= (dF[0] * F[0] + dF[2] * F[1]);
u0[1] -= (dF[1] * F[0] + dF[3] * F[1]);
u0[1] -= (dF[1] * F[0] + dF[3] * F[1]);*/
bool inbounds = true;
if (u0[0] < bounds[0])
......@@ -454,7 +454,7 @@ void FaceMesh::Smoothing()
continue;
}
Array<OneD, NekDouble> FN(2, 0.0);
/*Array<OneD, NekDouble> FN(2, 0.0);
pu0 = m_cadsurf->P(u0);
di = m_mesh->m_octree->Query(pu0);
for (int i = 0; i < nodesystem.size(); i++)
......@@ -477,7 +477,7 @@ void FaceMesh::Smoothing()
if (F[0] * F[0] + F[1] * F[1] < FN[0] * FN[0] + FN[1] * FN[1])
{
continue;
}
}*/
Array<OneD, NekDouble> l2 = m_cadsurf->P(u0);
(*nit)->Move(l2, m_id, u0);
......
......@@ -313,7 +313,7 @@ void InputMCF::Process()
ModuleKey(eProcessModule, "surfacemesh"), m_mesh));
////**** VolumeMesh ****////
/*mods.push_back(GetModuleFactory().CreateInstance(
mods.push_back(GetModuleFactory().CreateInstance(
ModuleKey(eProcessModule, "volumemesh"), m_mesh));
if(m_makeBL)
{
......@@ -321,11 +321,11 @@ void InputMCF::Process()
mods.back()->RegisterConfig("blthick",m_blthick);
mods.back()->RegisterConfig("bllayers",m_bllayers);
mods.back()->RegisterConfig("blprog",m_blprog);
}*/
}
}
////**** HOSurface ****////
/*mods.push_back(GetModuleFactory().CreateInstance(
mods.push_back(GetModuleFactory().CreateInstance(
ModuleKey(eProcessModule, "hosurface"), m_mesh));
if (m_surfopti)
{
......@@ -348,7 +348,7 @@ void InputMCF::Process()
}
////**** SPLIT BL ****////
/*if(m_splitBL)
if(m_splitBL)
{
mods.push_back(GetModuleFactory().CreateInstance(
ModuleKey(eProcessModule, "bl"), m_mesh));
......@@ -356,15 +356,13 @@ void InputMCF::Process()
mods.back()->RegisterConfig("surf",m_blsurfs);
mods.back()->RegisterConfig("nq",boost::lexical_cast<string>(m_mesh->m_nummode));
mods.back()->RegisterConfig("r",m_blprog);
}*/
}
for(int i = 0; i < mods.size(); i++)
{
mods[i]->SetDefaults();
mods[i]->Process();
}
m_mesh->m_expDim = 2;
}
}
}
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