43 #include <boost/algorithm/string.hpp> 44 #include <boost/filesystem.hpp> 50 namespace NekMeshUtils
54 "oce", CADSystemOCE::create,
"Uses OCE as cad engine");
62 if (modShape.IsBound(S))
70 bool CADSystemOCE::LoadCAD()
72 Handle(XSControl_WorkSession) WS;
73 Handle(Interface_InterfaceModel) Model;
74 Handle(XSControl_TransferReader) TR;
75 Handle(Transfer_TransientProcess) TP;
76 Handle(XCAFDoc_ShapeTool) STool;
78 bool fromStep =
false;
80 if (m_naca.size() == 0)
84 string ext = boost::filesystem::extension(m_name);
86 if (boost::iequals(ext,
".geo"))
88 shape = BuildGeo(m_name);
93 STEPCAFControl_Reader readerCAF;
94 readerCAF.SetNameMode(
true);
95 readerCAF.SetLayerMode(
true);
96 readerCAF.SetColorMode(
true);
98 Handle(TDocStd_Document) document =
99 new TDocStd_Document(Storage::Version());
100 readerCAF.ReadFile(m_name.c_str());
101 readerCAF.Transfer(document);
103 STEPControl_Reader reader = readerCAF.Reader();
107 TR = WS->TransferReader();
108 TP = TR->TransientProcess();
109 XCAFDoc_DocumentTool::ShapeTool(document->Main());
111 reader.NbRootsForTransfer();
112 reader.TransferRoots();
113 shape = reader.OneShape();
123 shape = BuildNACA(m_name);
126 TopExp_Explorer explr;
128 TopTools_DataMapOfShapeShape modShape;
132 BRepBuilderAPI_Sewing sew(1e-1);
134 for (explr.Init(shape, TopAbs_FACE); explr.More(); explr.Next())
136 sew.Add(explr.Current());
141 for (explr.Init(shape, TopAbs_FACE); explr.More(); explr.Next())
143 if (sew.IsModified(explr.Current()))
145 modShape.Bind(explr.Current(), sew.Modified(explr.Current()));
149 shape = sew.SewedShape();
152 for (explr.Init(shape, TopAbs_SHELL); explr.More(); explr.Next())
158 "Was not able to form a topological water tight shell");
162 for (explr.Init(shape, TopAbs_VERTEX); explr.More(); explr.Next())
164 TopoDS_Shape v = explr.Current();
165 if (mapOfVerts.Contains(v))
169 int i = mapOfVerts.Add(v);
176 for (explr.Init(shape, TopAbs_EDGE); explr.More(); explr.Next())
178 TopoDS_Shape e = explr.Current().Oriented(TopAbs_FORWARD);
179 if (mapOfEdges.Contains(e))
184 if (!BRep_Tool::Degenerated(TopoDS::Edge(e)))
186 int i = mapOfEdges.Add(e);
191 for (explr.Init(shape, TopAbs_FACE); explr.More(); explr.Next())
193 TopoDS_Shape f = explr.Current();
194 ASSERTL0(!mapOfFaces.Contains(f),
"duplicated faces");
195 int i = mapOfFaces.Add(f);
203 int nb = Model->NbEntities();
204 for (
int i = 1; i <= nb; i++)
206 if (!Model->Value(i)->DynamicType()->SubType(
207 "StepRepr_RepresentationItem"))
210 Handle(StepRepr_RepresentationItem) enti =
211 Handle(StepRepr_RepresentationItem)::DownCast(Model->Value(i));
212 Handle(TCollection_HAsciiString)
name = enti->Name();
217 Handle(Transfer_Binder) binder = TP->Find(Model->Value(i));
218 if (binder.IsNull() || !binder->HasResult())
221 TopoDS_Shape S = TransferBRep::ShapeResult(TP, binder);
226 if (S.ShapeType() == TopAbs_FACE)
228 string s(name->ToCString());
230 if (mapOfFaces.Contains(S))
232 int id = mapOfFaces.FindIndex(S);
234 m_surfs[id]->SetName(s);
239 if (mapOfFaces.Contains(S))
241 int id = mapOfFaces.FindIndex(S);
242 m_surfs[id]->SetName(s);
254 for (
int i = 1; i <= mapOfFaces.Extent(); i++)
256 TopoDS_Shape face = mapOfFaces.FindKey(i).Oriented(TopAbs_FORWARD);
258 TopTools_IndexedMapOfShape localEdges;
259 TopExp::MapShapes(face, TopAbs_EDGE, localEdges);
261 for (
int j = 1; j <= localEdges.Extent(); j++)
263 TopoDS_Shape edge = localEdges.FindKey(j);
264 if (BRep_Tool::Degenerated(TopoDS::Edge(edge)))
268 BRep_Tool::UVPoints(TopoDS::Edge(edge), TopoDS::Face(face), p1,
271 m_verts[mapOfVerts.FindIndex(TopExp::FirstVertex(
272 TopoDS::Edge(edge), Standard_True))]
273 ->SetDegen(i, m_surfs[i], (p1.X() + p2.X()) / 2.0,
274 (p1.Y() + p2.Y()) / 2.0);
282 for (
auto &i : m_curves)
284 ASSERTL0(i.second->GetAdjSurf().size() == 2,
285 "topolgy error found, surface not closed");
297 void CADSystemOCE::AddVert(
int i, TopoDS_Shape in)
301 std::static_pointer_cast<
CADVertOCE>(newVert)->Initialise(i, in);
303 m_verts[i] = newVert;
306 void CADSystemOCE::AddCurve(
int i, TopoDS_Shape in)
309 std::static_pointer_cast<
CADCurveOCE>(newCurve)->Initialise(i, in);
311 TopoDS_Vertex fv = TopExp::FirstVertex(TopoDS::Edge(in));
312 TopoDS_Vertex lv = TopExp::LastVertex(TopoDS::Edge(in));
314 vector<CADVertSharedPtr> vs;
315 vs.push_back(m_verts[mapOfVerts.FindIndex(fv)]);
316 vs.push_back(m_verts[mapOfVerts.FindIndex(lv)]);
317 newCurve->SetVert(vs);
319 m_curves[i] = newCurve;
322 void CADSystemOCE::AddSurf(
int i, TopoDS_Shape in)
325 std::static_pointer_cast<
CADSurfOCE>(newSurf)->Initialise(i, in);
328 TopoDS_Shape face = in.Oriented(TopAbs_FORWARD);
329 TopTools_IndexedMapOfShape mapOfWires;
330 TopExp::MapShapes(face, TopAbs_WIRE, mapOfWires);
331 vector<EdgeLoopSharedPtr> edgeloops;
333 for (
int j = 1; j <= mapOfWires.Extent(); j++)
337 TopoDS_Shape wire = mapOfWires.FindKey(j);
339 BRepTools_WireExplorer exp;
341 exp.Init(TopoDS::Wire(wire));
345 TopoDS_Shape edge = exp.Current();
347 if (mapOfEdges.Contains(edge))
349 int e = mapOfEdges.FindIndex(edge);
350 edgeloop->edges.push_back(m_curves[e]);
351 edgeloop->edgeo.push_back(exp.Orientation() == TopAbs_FORWARD
357 edgeloops.push_back(edgeloop);
361 for (
int k = 0; k < edgeloops.size(); k++)
363 tote += edgeloops[k]->edges.size();
366 ASSERTL0(tote != 1,
"cannot handle periodic curves");
368 CADSurf::OrientateEdges(newSurf, edgeloops);
369 newSurf->SetEdges(edgeloops);
372 for (
int k = 0; k < edgeloops.size(); k++)
374 for (
int j = 0; j < edgeloops[k]->edges.size(); j++)
376 edgeloops[k]->edges[j]->SetAdjSurf(
377 make_pair(newSurf, edgeloops[k]->edgeo[j]));
381 m_surfs[i] = newSurf;
387 bound[0] = numeric_limits<double>::max();
388 bound[1] = numeric_limits<double>::min();
389 bound[2] = numeric_limits<double>::max();
390 bound[3] = numeric_limits<double>::min();
391 bound[4] = numeric_limits<double>::max();
392 bound[5] = numeric_limits<double>::min();
394 for (
int i = 1; i <= m_curves.size(); i++)
399 bound[0] = min(bound[0], min(ends[0], ends[3]));
400 bound[1] = max(bound[1], max(ends[0], ends[3]));
402 bound[2] = min(bound[2], min(ends[1], ends[4]));
403 bound[3] = max(bound[3], max(ends[1], ends[4]));
405 bound[4] = min(bound[4], min(ends[2], ends[5]));
406 bound[5] = max(bound[5], max(ends[2], ends[5]));
412 TopoDS_Shape CADSystemOCE::BuildNACA(
string naca)
414 ASSERTL0(naca.length() == 4,
"not a 4 digit code");
415 vector<NekDouble> data;
416 ParseUtils::GenerateVector(m_naca, data);
417 ASSERTL0(data.size() == 5,
"not a vaild domain");
419 int n = std::stoi(naca);
430 for (
int i = 0; i < np; i++)
432 xc[i] = (1.0 - cos(i * dtheta)) / 2.0;
436 for (
int i = 0; i < np; i++)
440 yc[i] = M / P / P * (2.0 * P * xc[i] - xc[i] * xc[i]);
441 dyc[i] = 2.0 * M / P / P * (P - xc[i]);
445 yc[i] = M / (1.0 -
P) / (1.0 - P) *
446 (1.0 - 2.0 * P + 2.0 * P * xc[i] - xc[i] * xc[i]);
447 dyc[i] = 2.0 * M / (1.0 -
P) / (1.0 - P) * (P - xc[i]);
452 for (
int i = 0; i < np; i++)
455 T / 0.2 * (0.2969 * sqrt(xc[i]) - 0.1260 * xc[i] -
456 0.3516 * xc[i] * xc[i] + 0.2843 * xc[i] * xc[i] * xc[i] -
457 0.1015 * xc[i] * xc[i] * xc[i] * xc[i]);
462 for (
int i = np - 1; i >= 0; i--, l++)
466 x[l] = xc[i] - yt[i] * sin(theta);
467 y[l] = yc[i] + yt[i] * cos(theta);
469 for (
int i = 1; i < np; i++)
473 x[i + np - 1] = xc[i] + yt[i] * sin(theta);
474 y[i + np - 1] = yc[i] - yt[i] * cos(theta);
477 TColgp_Array1OfPnt pointArray(0, 2 * np - 2);
479 for (
int i = 0; i < 2 * np - 1; i++)
481 pointArray.SetValue(i, gp_Pnt(x[i] * 1000.0, y[i] * 1000.0, 0.0));
484 GeomAPI_PointsToBSpline spline(pointArray);
485 Handle(Geom_BSplineCurve) curve = spline.Curve();
487 BRepBuilderAPI_MakeEdge areoEdgeBuilder(curve);
488 TopoDS_Edge aeroEdge = areoEdgeBuilder.Edge();
489 BRepBuilderAPI_MakeEdge aeroTEBuilder(
490 gp_Pnt(x[0] * 1000.0, y[0] * 1000.0, 0.0),
491 gp_Pnt(x[2 * np - 2] * 1000.0, y[2 * np - 2] * 1000.0, 0.0));
492 TopoDS_Edge TeEdge = aeroTEBuilder.Edge();
494 BRepBuilderAPI_MakeWire aeroWireBuilder(aeroEdge, TeEdge);
495 TopoDS_Wire aeroWire = aeroWireBuilder.Wire();
498 gp_Ax1 rotAx(gp_Pnt(500.0, 0.0, 0.0), gp_Dir(gp_Vec(0.0, 0.0, -1.0)));
499 transform.SetRotation(rotAx, data[4] / 180.0 * M_PI);
500 TopLoc_Location mv(transform);
503 BRepBuilderAPI_MakeEdge domInlBuilder(
504 gp_Pnt(data[0] * 1000.0, data[1] * 1000.0, 0.0),
505 gp_Pnt(data[0] * 1000.0, data[3] * 1000.0, 0.0));
506 TopoDS_Edge inlEdge = domInlBuilder.Edge();
508 BRepBuilderAPI_MakeEdge domTopBuilder(
509 gp_Pnt(data[0] * 1000.0, data[3] * 1000.0, 0.0),
510 gp_Pnt(data[2] * 1000.0, data[3] * 1000.0, 0.0));
511 TopoDS_Edge topEdge = domTopBuilder.Edge();
513 BRepBuilderAPI_MakeEdge domOutBuilder(
514 gp_Pnt(data[2] * 1000.0, data[3] * 1000.0, 0.0),
515 gp_Pnt(data[2] * 1000.0, data[1] * 1000.0, 0.0));
516 TopoDS_Edge outEdge = domOutBuilder.Edge();
518 BRepBuilderAPI_MakeEdge domBotBuilder(
519 gp_Pnt(data[2] * 1000.0, data[1] * 1000.0, 0.0),
520 gp_Pnt(data[0] * 1000.0, data[1] * 1000.0, 0.0));
521 TopoDS_Edge botEdge = domBotBuilder.Edge();
523 BRepBuilderAPI_MakeWire domWireBuilder(inlEdge, topEdge, outEdge, botEdge);
524 TopoDS_Wire domWire = domWireBuilder.Wire();
526 BRepBuilderAPI_MakeFace face(domWire,
true);
529 ShapeFix_Face sf(face.Face());
535 TopoDS_Shape CADSystemOCE::BuildGeo(
string geo)
540 map<int, string> points;
541 map<int, string> lines;
542 map<int, string> splines;
543 map<int, string> bsplines;
544 map<int, string> circles;
545 map<int, string> ellipses;
546 map<int, string> loops;
547 map<int, string> surfs;
556 boost::erase_all(fline,
"\r");
558 vector<string> tmp1, tmp2;
559 boost::split(tmp1, fline, boost::is_any_of(
"//"));
562 if (!boost::contains(fline,
";"))
568 fline = flinetmp + fline;
571 boost::split(tmp1, fline, boost::is_any_of(
"="));
573 boost::split(tmp2, tmp1[0], boost::is_any_of(
"("));
575 string type = tmp2[0];
576 boost::erase_all(tmp2[1],
")");
577 boost::erase_all(tmp2[1],
" ");
578 int id = std::stoi(tmp2[1]);
580 boost::erase_all(tmp1[1],
" ");
581 boost::erase_all(tmp1[1],
"{");
582 boost::erase_all(tmp1[1],
"}");
583 boost::erase_all(tmp1[1],
";");
585 string var = tmp1[1];
587 if (boost::iequals(type,
"Point"))
591 else if (boost::iequals(type,
"Line"))
595 else if (boost::iequals(type,
"Spline"))
599 else if (boost::iequals(type,
"BSpline"))
603 else if (boost::iequals(type,
"Circle"))
607 else if (boost::iequals(type,
"Ellipse"))
611 else if (boost::iequals(type,
"Line Loop"))
616 boost::erase_all(var,
"-");
619 else if (boost::iequals(type,
"Plane Surface"))
625 cout <<
"not sure " << type << endl;
629 map<int, string>::iterator it;
632 map<int, gp_Pnt> cPoints;
633 for (it = points.begin(); it != points.end(); it++)
635 vector<NekDouble> data;
636 ParseUtils::GenerateVector(it->second, data);
639 gp_Pnt(data[0] * 1000.0, data[1] * 1000.0, data[2] * 1000.0);
643 map<int, TopoDS_Edge> cEdges;
644 for (it = lines.begin(); it != lines.end(); it++)
646 vector<unsigned int> data;
647 ParseUtils::GenerateVector(it->second, data);
648 BRepBuilderAPI_MakeEdge em(cPoints[data[0]], cPoints[data[1]]);
649 cEdges[it->first] = em.Edge();
651 for (it = splines.begin(); it != splines.end(); it++)
653 vector<unsigned int> data;
654 ParseUtils::GenerateVector(it->second, data);
656 TColgp_Array1OfPnt pointArray(0, data.size() - 1);
658 for (
int i = 0; i < data.size(); i++)
660 pointArray.SetValue(i, cPoints[data[i]]);
662 GeomAPI_PointsToBSpline spline(pointArray);
663 Handle(Geom_BSplineCurve) curve = spline.Curve();
665 BRepBuilderAPI_MakeEdge em(curve);
666 cEdges[it->first] = em.Edge();
668 for (it = bsplines.begin(); it != bsplines.end(); it++)
670 vector<unsigned int> data;
671 ParseUtils::GenerateVector(it->second, data);
673 TColgp_Array1OfPnt pointArray(0, data.size() - 1);
675 for (
int i = 0; i < data.size(); i++)
677 pointArray.SetValue(i, cPoints[data[i]]);
679 Handle(Geom_BezierCurve) curve =
new Geom_BezierCurve(pointArray);
681 BRepBuilderAPI_MakeEdge em(curve);
682 cEdges[it->first] = em.Edge();
684 for (it = circles.begin(); it != circles.end(); it++)
686 vector<unsigned int> data;
687 ParseUtils::GenerateVector(it->second, data);
689 ASSERTL0(data.size() == 3,
"Wrong definition of circle arc");
690 gp_Pnt start = cPoints[data[0]];
691 gp_Pnt centre = cPoints[data[1]];
692 gp_Pnt end = cPoints[data[2]];
696 ASSERTL0(fabs(r1 - r2) < 1e-7,
"Non-matching radii");
699 c.SetLocation(centre);
701 Handle(Geom_Circle) gc =
new Geom_Circle(c);
703 ShapeAnalysis_Curve sac;
705 sac.Project(gc, start, 1e-8, start, p1);
706 sac.Project(gc, end, 1e-8, end, p2);
709 if ((p1 > p2) ^ (fabs(p2 - p1) > M_PI))
714 Handle(Geom_TrimmedCurve) tc =
new Geom_TrimmedCurve(gc, p1, p2,
false);
716 BRepBuilderAPI_MakeEdge em(tc);
718 cEdges[it->first] = em.Edge();
720 for (it = ellipses.begin(); it != ellipses.end(); it++)
722 vector<unsigned int> data;
723 ParseUtils::GenerateVector(it->second, data);
725 ASSERTL0(data.size() == 4,
"Wrong definition of ellipse arc");
726 gp_Pnt start = cPoints[data[0]];
727 gp_Pnt centre = cPoints[data[1]];
729 gp_Pnt end = cPoints[data[3]];
731 NekDouble major = start.Distance(centre);
733 gp_Vec v1(centre, start);
734 gp_Vec v2(centre, end);
736 gp_Vec vx(1.0, 0.0, 0.0);
744 v2.Rotate(gp_Ax1(), angle);
745 NekDouble minor = fabs(v2.Y() / sqrt(1.0 - v2.X() * v2.X() / (major * major)));
748 e.SetLocation(centre);
749 e.SetMajorRadius(major);
750 e.SetMinorRadius(minor);
751 e.Rotate(e.Axis(), angle);
752 Handle(Geom_Ellipse) ge =
new Geom_Ellipse(e);
754 ShapeAnalysis_Curve sac;
756 sac.Project(ge, start, 1e-8, start, p1);
757 sac.Project(ge, end, 1e-8, end, p2);
760 if (fabs(p2 - p1) > M_PI)
765 Handle(Geom_TrimmedCurve) tc =
new Geom_TrimmedCurve(ge, p1, p2,
true);
767 BRepBuilderAPI_MakeEdge em(tc);
769 cEdges[it->first] = em.Edge();
773 map<int, TopoDS_Wire> cWires;
774 for (it = loops.begin(); it != loops.end(); it++)
776 vector<unsigned int> data;
777 ParseUtils::GenerateVector(it->second, data);
778 BRepBuilderAPI_MakeWire wm;
779 for (
int i = 0; i < data.size(); i++)
781 wm.Add(cEdges[data[i]]);
783 cWires[it->first] = wm.Wire();
788 ASSERTL0(surfs.size() == 1,
"more than 1 surf");
790 vector<unsigned int> data;
791 ParseUtils::GenerateVector(it->second, data);
792 BRepBuilderAPI_MakeFace face(cWires[data[0]],
true);
793 for (
int i = 1; i < data.size(); i++)
795 face.Add(cWires[data[i]]);
798 ASSERTL0(face.Error() == BRepBuilderAPI_FaceDone,
"build geo failed");
800 ShapeFix_Face sf(face.Face());
CADCurveFactory & GetCADCurveFactory()
std::shared_ptr< CADSurf > CADSurfSharedPtr
#define ASSERTL0(condition, msg)
std::shared_ptr< EdgeLoop > EdgeLoopSharedPtr
std::shared_ptr< CADVert > CADVertSharedPtr
struct which descibes a collection of cad edges which are a loop on the cad surface ...
CADVertFactory & GetCADVertFactory()
void filterModShape(TopTools_DataMapOfShapeShape &modShape, TopoDS_Shape &S)
CADSurfFactory & GetCADSurfFactory()
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
std::shared_ptr< CADCurve > CADCurveSharedPtr
EngineFactory & GetEngineFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.