Commit 41981264 authored by Michael Turner's avatar Michael Turner

working

parent 4a623001
......@@ -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()
......
......@@ -72,10 +72,10 @@ bool CADSystemOCE::LoadCAD()
cout << "assuming " << m_name << " is a 4 digit naca code" << endl;
shape = BuildNACA(m_name);
STEPControl_Writer writer;
/*STEPControl_Writer writer;
writer.Transfer(shape,STEPControl_ShellBasedSurfaceModel);
writer.Write("test.stp");
exit(-1);
exit(-1);*/
}
// faces and verts can be extracted straight from shape
......@@ -375,15 +375,17 @@ TopoDS_Shape CADSystemOCE::BuildNACA(string naca)
n/=10;
NekDouble M = (n%10)/100.0;
Array<OneD, NekDouble> xc(100);
NekDouble dtheta = M_PI/99;
for(int i = 0; i < 100; i++)
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(100), dyc(100);
for(int i = 0; i < 100; i++)
Array<OneD, NekDouble> yc(np), dyc(np);
for(int i = 0; i < np; i++)
{
if(xc[i] < P)
{
......@@ -397,8 +399,8 @@ TopoDS_Shape CADSystemOCE::BuildNACA(string naca)
}
}
Array<OneD, NekDouble> yt(100);
for(int i = 0; i < 100; 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]
......@@ -407,9 +409,9 @@ TopoDS_Shape CADSystemOCE::BuildNACA(string naca)
-0.1015 * xc[i] * xc[i] * xc[i] * xc[i]);
}
Array<OneD, NekDouble> x(199), y(199);
Array<OneD, NekDouble> x(2*np-1), y(2*np-1);
int l = 0;
for(int i = 99; i >= 0; i--, l++)
for(int i = np-1; i >= 0; i--, l++)
{
NekDouble theta = atan(dyc[i]);
......@@ -417,41 +419,39 @@ TopoDS_Shape CADSystemOCE::BuildNACA(string naca)
y[l] = yc[i] + yt[i] * cos(theta);
}
for(int i = 1; i < 100; i++)
for(int i = 1; i < np; i++)
{
NekDouble theta = atan(dyc[i]);
x[i+99] = xc[i] + yt[i] * sin(theta);
y[i+99] = yc[i] - yt[i] * cos(theta);
x[i+np-1] = xc[i] + yt[i] * sin(theta);
y[i+np-1] = yc[i] - yt[i] * cos(theta);
}
Handle(TColgp_HArray1OfPnt) pointArray = new TColgp_HArray1OfPnt(0,198);
TColgp_Array1OfPnt pointArray(0,2*np-2);
for(int i = 0; i < 199; i++)
for(int i = 0; i < 2*np-1; i++)
{
pointArray->SetValue(i,gp_Pnt(x[i],y[i],0.0));
pointArray.SetValue(i,gp_Pnt(x[i]*1000.0,y[i]*1000.0,0.0));
}
TopLoc_Location loc;
GeomAPI_Interpolate spline(pointArray,false,1e-8);
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],y[0],0.0), gp_Pnt(x[198],y[198],0.0));
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(-2.0,-2.0,0.0), gp_Pnt(-2.0,2.0,0.0));
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(-2.0,2.0,0.0), gp_Pnt(5.0,2.0,0.0));
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(5.0,2.0,0.0), gp_Pnt(5.0,-2.0,0.0));
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(5.0,-2.0,0.0), gp_Pnt(-2.0,-2.0,0.0));
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);
......
......@@ -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,8 +71,9 @@
#include <ShapeAnalysis_Curve.hxx>
#include <Standard_Macro.hxx>
#include <GeomAPI_Interpolate.hxx>
#include <TColgp_HArray1OfPnt.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>
......
......@@ -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