47 namespace NekMeshUtils
51 "oce", CADSystemOCE::create,
"Uses OCE as cad engine");
53 bool CADSystemOCE::LoadCAD()
55 if (m_naca.size() == 0)
59 STEPControl_Reader reader;
60 reader = STEPControl_Reader();
61 reader.ReadFile(m_name.c_str());
62 reader.NbRootsForTransfer();
63 reader.TransferRoots();
64 shape = reader.OneShape();
72 shape = BuildNACA(m_name);
75 TopExp_Explorer explr;
78 for (explr.Init(shape, TopAbs_VERTEX); explr.More(); explr.Next())
80 TopoDS_Shape v = explr.Current();
81 if (mapOfVerts.Contains(v))
85 int i = mapOfVerts.Add(v);
92 for (explr.Init(shape, TopAbs_EDGE); explr.More(); explr.Next())
94 TopoDS_Shape e = explr.Current().Oriented(TopAbs_FORWARD);
95 if (mapOfEdges.Contains(e))
99 BRepAdaptor_Curve curve = BRepAdaptor_Curve(TopoDS::Edge(e));
100 if (curve.GetType() != 7)
102 int i = mapOfEdges.Add(e);
107 for (explr.Init(shape, TopAbs_FACE); explr.More(); explr.Next())
109 TopoDS_Shape f = explr.Current();
110 ASSERTL0(!mapOfFaces.Contains(f),
"duplicated faces");
111 int i = mapOfFaces.Add(f);
117 for (
int i = 1; i <= mapOfFaces.Extent(); i++)
119 TopoDS_Shape face = mapOfFaces.FindKey(i).Oriented(TopAbs_FORWARD);
121 TopTools_IndexedMapOfShape localEdges;
122 TopExp::MapShapes(face, TopAbs_EDGE, localEdges);
124 for (
int j = 1; j <= localEdges.Extent(); j++)
126 TopoDS_Shape edge = localEdges.FindKey(j);
127 if (BRep_Tool::Degenerated(TopoDS::Edge(edge)))
131 BRep_Tool::UVPoints(TopoDS::Edge(edge), TopoDS::Face(face), p1,
134 m_verts[mapOfVerts.FindIndex(TopExp::FirstVertex(
135 TopoDS::Edge(edge), Standard_True))]
136 ->SetDegen(i, m_surfs[i], (p1.X() + p2.X()) / 2.0,
137 (p1.Y() + p2.Y()) / 2.0);
146 for (it = m_curves.begin(); it != m_curves.end(); it++)
148 ASSERTL0(it->second->GetAdjSurf().size() == 2,
149 "curve is not joined to 2 surfaces");
156 void CADSystemOCE::AddVert(
int i, TopoDS_Shape in)
160 boost::static_pointer_cast<
CADVertOCE>(newVert)->Initialise(i, in);
162 m_verts[i] = newVert;
165 void CADSystemOCE::AddCurve(
int i, TopoDS_Shape in)
168 boost::static_pointer_cast<
CADCurveOCE>(newCurve)->Initialise(i, in);
170 TopoDS_Vertex fv = TopExp::FirstVertex(TopoDS::Edge(in));
171 TopoDS_Vertex lv = TopExp::LastVertex(TopoDS::Edge(in));
173 vector<CADVertSharedPtr> vs;
174 vs.push_back(m_verts[mapOfVerts.FindIndex(fv)]);
175 vs.push_back(m_verts[mapOfVerts.FindIndex(lv)]);
176 newCurve->SetVert(vs);
178 m_curves[i] = newCurve;
181 void CADSystemOCE::AddSurf(
int i, TopoDS_Shape in)
184 boost::static_pointer_cast<
CADSurfOCE>(newSurf)->Initialise(i, in);
187 TopoDS_Shape face = in.Oriented(TopAbs_FORWARD);
188 TopTools_IndexedMapOfShape mapOfWires;
189 TopExp::MapShapes(face, TopAbs_WIRE, mapOfWires);
190 vector<EdgeLoopSharedPtr> edgeloops;
192 for (
int j = 1; j <= mapOfWires.Extent(); j++)
196 TopoDS_Shape wire = mapOfWires.FindKey(j);
198 BRepTools_WireExplorer exp;
200 exp.Init(TopoDS::Wire(wire));
204 TopoDS_Shape edge = exp.Current();
206 if (mapOfEdges.Contains(edge))
208 int e = mapOfEdges.FindIndex(edge);
209 edgeloop->edges.push_back(m_curves[e]);
210 edgeloop->edgeo.push_back(exp.Orientation() == TopAbs_FORWARD
216 edgeloops.push_back(edgeloop);
220 for (
int k = 0; k < edgeloops.size(); k++)
222 tote += edgeloops[k]->edges.size();
225 ASSERTL0(tote != 1,
"cannot handle periodic curves");
227 CADSurf::OrientateEdges(newSurf, edgeloops);
228 newSurf->SetEdges(edgeloops);
231 for (
int k = 0; k < edgeloops.size(); k++)
233 for (
int j = 0; j < edgeloops[k]->edges.size(); j++)
235 edgeloops[k]->edges[j]->SetAdjSurf(
236 make_pair(newSurf, edgeloops[k]->edgeo[j]));
240 m_surfs[i] = newSurf;
246 bound[0] = numeric_limits<double>::max();
247 bound[1] = numeric_limits<double>::min();
248 bound[2] = numeric_limits<double>::max();
249 bound[3] = numeric_limits<double>::min();
250 bound[4] = numeric_limits<double>::max();
251 bound[5] = numeric_limits<double>::min();
253 for (
int i = 1; i <= m_curves.size(); i++)
258 bound[0] = min(bound[0], min(ends[0], ends[3]));
259 bound[1] = max(bound[1], max(ends[0], ends[3]));
261 bound[2] = min(bound[2], min(ends[1], ends[4]));
262 bound[3] = max(bound[3], max(ends[1], ends[4]));
264 bound[4] = min(bound[4], min(ends[2], ends[5]));
265 bound[5] = max(bound[5], max(ends[2], ends[5]));
271 TopoDS_Shape CADSystemOCE::BuildNACA(
string naca)
273 ASSERTL0(naca.length() == 4,
"not a 4 digit code");
274 vector<NekDouble> data;
275 ParseUtils::GenerateUnOrderedVector(m_naca.c_str(), data);
276 ASSERTL0(data.size() == 5,
"not a vaild domain");
278 int n = boost::lexical_cast<
int>(naca);
289 for (
int i = 0; i < np; i++)
291 xc[i] = (1.0 - cos(i * dtheta)) / 2.0;
295 for (
int i = 0; i < np; i++)
299 yc[i] = M / P / P * (2.0 * P * xc[i] - xc[i] * xc[i]);
300 dyc[i] = 2.0 * M / P / P * (P - xc[i]);
304 yc[i] = M / (1.0 -
P) / (1.0 - P) *
305 (1.0 - 2.0 * P + 2.0 * P * xc[i] - xc[i] * xc[i]);
306 dyc[i] = 2.0 * M / (1.0 -
P) / (1.0 - P) * (P - xc[i]);
311 for (
int i = 0; i < np; i++)
314 T / 0.2 * (0.2969 * sqrt(xc[i]) - 0.1260 * xc[i] -
315 0.3516 * xc[i] * xc[i] + 0.2843 * xc[i] * xc[i] * xc[i] -
316 0.1015 * xc[i] * xc[i] * xc[i] * xc[i]);
321 for (
int i = np - 1; i >= 0; i--, l++)
325 x[l] = xc[i] - yt[i] * sin(theta);
326 y[l] = yc[i] + yt[i] * cos(theta);
328 for (
int i = 1; i < np; i++)
332 x[i + np - 1] = xc[i] + yt[i] * sin(theta);
333 y[i + np - 1] = yc[i] - yt[i] * cos(theta);
336 TColgp_Array1OfPnt pointArray(0, 2 * np - 2);
338 for (
int i = 0; i < 2 * np - 1; i++)
340 pointArray.SetValue(i, gp_Pnt(x[i] * 1000.0, y[i] * 1000.0, 0.0));
343 GeomAPI_PointsToBSpline spline(pointArray);
344 Handle(Geom_BSplineCurve) curve = spline.Curve();
346 BRepBuilderAPI_MakeEdge areoEdgeBuilder(curve);
347 TopoDS_Edge aeroEdge = areoEdgeBuilder.Edge();
348 BRepBuilderAPI_MakeEdge aeroTEBuilder(
349 gp_Pnt(x[0] * 1000.0, y[0] * 1000.0, 0.0),
350 gp_Pnt(x[2 * np - 2] * 1000.0, y[2 * np - 2] * 1000.0, 0.0));
351 TopoDS_Edge TeEdge = aeroTEBuilder.Edge();
353 BRepBuilderAPI_MakeWire aeroWireBuilder(aeroEdge, TeEdge);
354 TopoDS_Wire aeroWire = aeroWireBuilder.Wire();
357 gp_Ax1 rotAx(gp_Pnt(500.0, 0.0, 0.0), gp_Dir(gp_Vec(0.0, 0.0, -1.0)));
358 transform.SetRotation(rotAx, data[4] / 180.0 * M_PI);
359 TopLoc_Location mv(transform);
362 BRepBuilderAPI_MakeEdge domInlBuilder(
363 gp_Pnt(data[0] * 1000.0, data[1] * 1000.0, 0.0),
364 gp_Pnt(data[0] * 1000.0, data[3] * 1000.0, 0.0));
365 TopoDS_Edge inlEdge = domInlBuilder.Edge();
367 BRepBuilderAPI_MakeEdge domTopBuilder(
368 gp_Pnt(data[0] * 1000.0, data[3] * 1000.0, 0.0),
369 gp_Pnt(data[2] * 1000.0, data[3] * 1000.0, 0.0));
370 TopoDS_Edge topEdge = domTopBuilder.Edge();
372 BRepBuilderAPI_MakeEdge domOutBuilder(
373 gp_Pnt(data[2] * 1000.0, data[3] * 1000.0, 0.0),
374 gp_Pnt(data[2] * 1000.0, data[1] * 1000.0, 0.0));
375 TopoDS_Edge outEdge = domOutBuilder.Edge();
377 BRepBuilderAPI_MakeEdge domBotBuilder(
378 gp_Pnt(data[2] * 1000.0, data[1] * 1000.0, 0.0),
379 gp_Pnt(data[0] * 1000.0, data[1] * 1000.0, 0.0));
380 TopoDS_Edge botEdge = domBotBuilder.Edge();
382 BRepBuilderAPI_MakeWire domWireBuilder(inlEdge, topEdge, outEdge, botEdge);
383 TopoDS_Wire domWire = domWireBuilder.Wire();
385 BRepBuilderAPI_MakeFace face(domWire,
true);
388 ShapeFix_Face sf(face.Face());
CADCurveFactory & GetCADCurveFactory()
#define ASSERTL0(condition, msg)
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
CADVertFactory & GetCADVertFactory()
CADSurfFactory & GetCADSurfFactory()
boost::shared_ptr< CADSurf > CADSurfSharedPtr
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
EngineFactory & GetEngineFactory()
struct which descibes a collection of cad edges which are a loop on the cad surface ...
boost::shared_ptr< EdgeLoop > EdgeLoopSharedPtr
boost::shared_ptr< CADVert > CADVertSharedPtr
boost::shared_ptr< CADCurve > CADCurveSharedPtr
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.