Commit 486cdf2d authored by Michael Turner's avatar Michael Turner

Merge branch 'feature/2D-generator-NACA' into feature/2D-generator

parents bf129784 41981264
......@@ -130,9 +130,10 @@ Array<OneD, NekDouble> CADCurveOCE::D2(NekDouble t)
NekDouble CADCurveOCE::Curvature(NekDouble t)
{
Array<OneD, NekDouble> der = D2(t);
GeomLProp_CLProps d(m_c,2,1e-8);
d.SetParameter(t);
return der[6] * der[6] + der[7] * der[7] + der[8] * der[8];
return d.Curvature() * 1000.0;
}
Array<OneD, NekDouble> CADCurveOCE::Bounds()
......
......@@ -53,16 +53,29 @@ bool CADSystemOCE::LoadCAD()
{
cout << "trying " << m_name << endl;
// Takes step file and makes OpenCascade shape
STEPControl_Reader reader;
reader = STEPControl_Reader();
reader.ReadFile(m_name.c_str());
reader.NbRootsForTransfer();
reader.TransferRoots();
shape = reader.OneShape();
if (shape.IsNull())
if (m_name.find('.') != std::string::npos)
{
return false;
// Takes step file and makes OpenCascade shape
STEPControl_Reader reader;
reader = STEPControl_Reader();
reader.ReadFile(m_name.c_str());
reader.NbRootsForTransfer();
reader.TransferRoots();
shape = reader.OneShape();
if (shape.IsNull())
{
return false;
}
}
else
{
cout << "assuming " << m_name << " is a 4 digit naca code" << endl;
shape = BuildNACA(m_name);
/*STEPControl_Writer writer;
writer.Transfer(shape,STEPControl_ShellBasedSurfaceModel);
writer.Write("test.stp");
exit(-1);*/
}
// faces and verts can be extracted straight from shape
......@@ -351,5 +364,104 @@ Array<OneD, NekDouble> CADSystemOCE::GetBoundingBox()
return bound;
}
TopoDS_Shape CADSystemOCE::BuildNACA(string naca)
{
ASSERTL0(naca.length() == 4, "not a 4 digit code");
int n = boost::lexical_cast<int>(naca);
NekDouble T = (n%100) / 100.0;
n/=100;
NekDouble P = (n%10)/10.0;
n/=10;
NekDouble M = (n%10)/100.0;
int np = 25;
Array<OneD, NekDouble> xc(np);
NekDouble dtheta = M_PI/(np-1);
for(int i = 0; i < np; i++)
{
xc[i] = (1.0 - cos(i*dtheta)) / 2.0;
}
Array<OneD, NekDouble> yc(np), dyc(np);
for(int i = 0; i < np; i++)
{
if(xc[i] < P)
{
yc[i] = M / P / P * (2.0 * P * xc[i] - xc[i] * xc[i]);
dyc[i] = 2.0 * M / P / P * (P - xc[i]);
}
else
{
yc[i] = M / (1.0 - P) / (1.0 - P) * (1.0 - 2.0 * P + 2.0 * P * xc[i] - xc[i] * xc[i]);
dyc[i] = 2.0 * M / (1.0 - P) / (1.0 - P) * (P - xc[i]);
}
}
Array<OneD, NekDouble> yt(np);
for(int i = 0; i < np; i++)
{
yt[i] = T / 0.2 * ( 0.2969 * sqrt(xc[i])
-0.1260 * xc[i]
-0.3516 * xc[i] * xc[i]
+0.2843 * xc[i] * xc[i] * xc[i]
-0.1015 * xc[i] * xc[i] * xc[i] * xc[i]);
}
Array<OneD, NekDouble> x(2*np-1), y(2*np-1);
int l = 0;
for(int i = np-1; i >= 0; i--, l++)
{
NekDouble theta = atan(dyc[i]);
x[l] = xc[i] - yt[i] * sin(theta);
y[l] = yc[i] + yt[i] * cos(theta);
}
for(int i = 1; i < np; i++)
{
NekDouble theta = atan(dyc[i]);
x[i+np-1] = xc[i] + yt[i] * sin(theta);
y[i+np-1] = yc[i] - yt[i] * cos(theta);
}
TColgp_Array1OfPnt pointArray(0,2*np-2);
for(int i = 0; i < 2*np-1; i++)
{
pointArray.SetValue(i,gp_Pnt(x[i]*1000.0,y[i]*1000.0,0.0));
}
GeomAPI_PointsToBSpline spline(pointArray);
Handle(Geom_BSplineCurve) curve = spline.Curve();
BRepBuilderAPI_MakeEdge areoEdgeBuilder(curve);
TopoDS_Edge aeroEdge = areoEdgeBuilder.Edge();
BRepBuilderAPI_MakeEdge aeroTEBuilder(gp_Pnt(x[0]*1000.0,y[0]*1000.0,0.0), gp_Pnt(x[2*np-2]*1000.0,y[2*np-2]*1000.0,0.0));
TopoDS_Edge TeEdge = aeroTEBuilder.Edge();
BRepBuilderAPI_MakeWire aeroWireBuilder(aeroEdge, TeEdge);
TopoDS_Wire aeroWire = aeroWireBuilder.Wire();
BRepBuilderAPI_MakeEdge domInlBuilder(gp_Pnt(-2000.0,-2000.0,0.0), gp_Pnt(-2000.0,2000.0,0.0));
TopoDS_Edge inlEdge = domInlBuilder.Edge();
BRepBuilderAPI_MakeEdge domTopBuilder(gp_Pnt(-2000.0,2000.0,0.0), gp_Pnt(5000.0,2000.0,0.0));
TopoDS_Edge topEdge = domTopBuilder.Edge();
BRepBuilderAPI_MakeEdge domOutBuilder(gp_Pnt(5000.0,2000.0,0.0), gp_Pnt(5000.0,-2000.0,0.0));
TopoDS_Edge outEdge = domOutBuilder.Edge();
BRepBuilderAPI_MakeEdge domBotBuilder(gp_Pnt(5000.0,-2000.0,0.0), gp_Pnt(-2000.0,-2000.0,0.0));
TopoDS_Edge botEdge = domBotBuilder.Edge();
BRepBuilderAPI_MakeWire domWireBuilder(inlEdge, topEdge, outEdge, botEdge);
TopoDS_Wire domWire = domWireBuilder.Wire();
BRepBuilderAPI_MakeFace face(domWire, true);
face.Add(aeroWire);
return face.Face();
}
}
}
......@@ -72,6 +72,8 @@ private:
void AddCurve(int i, TopoDS_Shape in, int fv, int lv);
/// Function to add surface to CADSystem::m_surfs.
void AddSurf(int i, TopoDS_Shape in, std::vector<EdgeLoop> ein);
TopoDS_Shape BuildNACA(std::string naca);
/// OCC master object
TopoDS_Shape shape;
};
......
......@@ -57,6 +57,7 @@
#include <GProp_GProps.hxx>
#include <BRepGProp.hxx>
#include <Geom_TrimmedCurve.hxx>
#include <GeomLProp_CLProps.hxx>
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <GCPnts_AbscissaPoint.hxx>
#include <TopAbs_State.hxx>
......@@ -70,4 +71,12 @@
#include <ShapeAnalysis_Curve.hxx>
#include <Standard_Macro.hxx>
#include <GeomAPI_PointsToBSpline.hxx>
#include <Geom_BSplineCurve.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRepBuilderAPI_MakeWire.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <STEPControl_Writer.hxx>
#endif
......@@ -873,10 +873,10 @@ void Octree::CompileSourcePointList()
CADCurveSharedPtr curve = m_mesh->m_cad->GetCurve(i);
Array<OneD, NekDouble> bds = curve->Bounds();
int samples = 100;
int dt = (bds[1] - bds[0]) / (samples + 1);
NekDouble dt = (bds[1] - bds[0]) / (samples + 1);
for (int j = 0; j < samples; j++)
{
NekDouble t = bds[0] + dt * j;
NekDouble t = bds[0] + dt * j;
NekDouble C = curve->Curvature(t);
Array<OneD, NekDouble> loc = curve->P(t);
......@@ -886,6 +886,8 @@ void Octree::CompileSourcePointList()
if (C != 0.0)
{
cout << 1.0 / C << endl;
NekDouble del = 2.0 * (1.0 / C) * sqrt(m_eps * (2.0 - m_eps));
if (del > m_maxDelta)
......
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