Nektar++
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CADSystemOCE.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: CADSystem.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: cad object methods.
32 //
33 ////////////////////////////////////////////////////////////////////////////////
34 
36 
42 
43 #include <boost/algorithm/string.hpp>
44 #include <boost/filesystem.hpp>
45 
46 using namespace std;
47 
48 namespace Nektar
49 {
50 namespace NekMeshUtils
51 {
52 
53 std::string CADSystemOCE::key = GetEngineFactory().RegisterCreatorFunction(
54  "oce", CADSystemOCE::create, "Uses OCE as cad engine");
55 
56 void filterModShape(TopTools_DataMapOfShapeShape &modShape, TopoDS_Shape &S)
57 {
58  bool repeat = true;
59  while (repeat)
60  {
61  repeat = false;
62  if (modShape.IsBound(S))
63  {
64  repeat = true;
65  S = modShape.Find(S);
66  }
67  }
68 }
69 
70 bool CADSystemOCE::LoadCAD()
71 {
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;
77 
78  bool fromStep = false;
79 
80  if (m_naca.size() == 0)
81  {
82  // not a naca profile behave normally
83  // could be a geo
84  string ext = boost::filesystem::extension(m_name);
85 
86  if (boost::iequals(ext, ".geo"))
87  {
88  shape = BuildGeo(m_name);
89  }
90  else
91  {
92  // Takes step file and makes OpenCascade shape
93  STEPCAFControl_Reader readerCAF;
94  readerCAF.SetNameMode(true);
95  readerCAF.SetLayerMode(true);
96  readerCAF.SetColorMode(true);
97 
98  Handle(TDocStd_Document) document =
99  new TDocStd_Document(Storage::Version());
100  readerCAF.ReadFile(m_name.c_str());
101  readerCAF.Transfer(document);
102 
103  STEPControl_Reader reader = readerCAF.Reader();
104 
105  WS = reader.WS();
106  Model = WS->Model();
107  TR = WS->TransferReader();
108  TP = TR->TransientProcess();
109  XCAFDoc_DocumentTool::ShapeTool(document->Main());
110 
111  reader.NbRootsForTransfer();
112  reader.TransferRoots();
113  shape = reader.OneShape();
114  if (shape.IsNull())
115  {
116  return false;
117  }
118  fromStep = true;
119  }
120  }
121  else
122  {
123  shape = BuildNACA(m_name);
124  }
125 
126  TopExp_Explorer explr;
127 
128  TopTools_DataMapOfShapeShape modShape;
129 
130  if(!m_2d)
131  {
132  BRepBuilderAPI_Sewing sew(1e-1);
133 
134  for (explr.Init(shape, TopAbs_FACE); explr.More(); explr.Next())
135  {
136  sew.Add(explr.Current());
137  }
138 
139  sew.Perform();
140 
141  for (explr.Init(shape, TopAbs_FACE); explr.More(); explr.Next())
142  {
143  if (sew.IsModified(explr.Current()))
144  {
145  modShape.Bind(explr.Current(), sew.Modified(explr.Current()));
146  }
147  }
148 
149  shape = sew.SewedShape();
150 
151  int shell = 0;
152  for (explr.Init(shape, TopAbs_SHELL); explr.More(); explr.Next())
153  {
154  shell++;
155  }
156 
157  ASSERTL0(shell == 1,
158  "Was not able to form a topological water tight shell");
159  }
160 
161  // build map of verticies
162  for (explr.Init(shape, TopAbs_VERTEX); explr.More(); explr.Next())
163  {
164  TopoDS_Shape v = explr.Current();
165  if (mapOfVerts.Contains(v))
166  {
167  continue;
168  }
169  int i = mapOfVerts.Add(v);
170  AddVert(i, v);
171  }
172 
173  // For each face of the geometry, get the local edges which bound it. If
174  // they are valid (their type != 7), then add them to an edge map. This
175  // filters out the dummy edges which OCC uses.
176  for (explr.Init(shape, TopAbs_EDGE); explr.More(); explr.Next())
177  {
178  TopoDS_Shape e = explr.Current().Oriented(TopAbs_FORWARD);
179  if (mapOfEdges.Contains(e))
180  {
181  continue;
182  }
183 
184  if (!BRep_Tool::Degenerated(TopoDS::Edge(e)))
185  {
186  int i = mapOfEdges.Add(e);
187  AddCurve(i, e);
188  }
189  }
190 
191  for (explr.Init(shape, TopAbs_FACE); explr.More(); explr.Next())
192  {
193  TopoDS_Shape f = explr.Current();
194  ASSERTL0(!mapOfFaces.Contains(f), "duplicated faces");
195  int i = mapOfFaces.Add(f);
196 
197  AddSurf(i, f);
198  }
199 
200  // attemps to extract patch names from STEP file
201  if(fromStep)
202  {
203  int nb = Model->NbEntities();
204  for (int i = 1; i <= nb; i++)
205  {
206  if (!Model->Value(i)->DynamicType()->SubType(
207  "StepRepr_RepresentationItem"))
208  continue;
209 
210  Handle(StepRepr_RepresentationItem) enti =
211  Handle(StepRepr_RepresentationItem)::DownCast(Model->Value(i));
212  Handle(TCollection_HAsciiString) name = enti->Name();
213 
214  if (name->IsEmpty())
215  continue;
216 
217  Handle(Transfer_Binder) binder = TP->Find(Model->Value(i));
218  if (binder.IsNull() || !binder->HasResult())
219  continue;
220 
221  TopoDS_Shape S = TransferBRep::ShapeResult(TP, binder);
222 
223  if (S.IsNull())
224  continue;
225 
226  if (S.ShapeType() == TopAbs_FACE)
227  {
228  string s(name->ToCString());
229 
230  if (mapOfFaces.Contains(S))
231  {
232  int id = mapOfFaces.FindIndex(S);
233 
234  m_surfs[id]->SetName(s);
235  }
236  else
237  {
238  filterModShape(modShape, S);
239  if (mapOfFaces.Contains(S))
240  {
241  int id = mapOfFaces.FindIndex(S);
242  m_surfs[id]->SetName(s);
243  }
244  else
245  {
246  ASSERTL0(false, "Name error");
247  }
248  }
249  }
250  }
251  }
252 
253  // attempts to identify properties of the vertex on the degen edge
254  for (int i = 1; i <= mapOfFaces.Extent(); i++)
255  {
256  TopoDS_Shape face = mapOfFaces.FindKey(i).Oriented(TopAbs_FORWARD);
257 
258  TopTools_IndexedMapOfShape localEdges;
259  TopExp::MapShapes(face, TopAbs_EDGE, localEdges);
260 
261  for (int j = 1; j <= localEdges.Extent(); j++)
262  {
263  TopoDS_Shape edge = localEdges.FindKey(j);
264  if (BRep_Tool::Degenerated(TopoDS::Edge(edge)))
265  {
266  gp_Pnt2d p1, p2;
267 
268  BRep_Tool::UVPoints(TopoDS::Edge(edge), TopoDS::Face(face), p1,
269  p2);
270 
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);
275  }
276  }
277  }
278 
279  // This checks that all edges are bound by two surfaces, sanity check.
280  if (!m_2d)
281  {
282  for (auto &i : m_curves)
283  {
284  ASSERTL0(i.second->GetAdjSurf().size() == 2,
285  "topolgy error found, surface not closed");
286  }
287  }
288 
289  if(m_verbose)
290  {
291  Report();
292  }
293 
294  return true;
295 }
296 
297 void CADSystemOCE::AddVert(int i, TopoDS_Shape in)
298 {
300 
301  std::static_pointer_cast<CADVertOCE>(newVert)->Initialise(i, in);
302 
303  m_verts[i] = newVert;
304 }
305 
306 void CADSystemOCE::AddCurve(int i, TopoDS_Shape in)
307 {
309  std::static_pointer_cast<CADCurveOCE>(newCurve)->Initialise(i, in);
310 
311  TopoDS_Vertex fv = TopExp::FirstVertex(TopoDS::Edge(in));
312  TopoDS_Vertex lv = TopExp::LastVertex(TopoDS::Edge(in));
313 
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);
318 
319  m_curves[i] = newCurve;
320 }
321 
322 void CADSystemOCE::AddSurf(int i, TopoDS_Shape in)
323 {
325  std::static_pointer_cast<CADSurfOCE>(newSurf)->Initialise(i, in);
326 
327  // do the exploration on forward oriented
328  TopoDS_Shape face = in.Oriented(TopAbs_FORWARD);
329  TopTools_IndexedMapOfShape mapOfWires;
330  TopExp::MapShapes(face, TopAbs_WIRE, mapOfWires);
331  vector<EdgeLoopSharedPtr> edgeloops;
332  // now we acutally analyse the loops for cad building
333  for (int j = 1; j <= mapOfWires.Extent(); j++)
334  {
336 
337  TopoDS_Shape wire = mapOfWires.FindKey(j);
338 
339  BRepTools_WireExplorer exp;
340 
341  exp.Init(TopoDS::Wire(wire));
342 
343  while (exp.More())
344  {
345  TopoDS_Shape edge = exp.Current();
346 
347  if (mapOfEdges.Contains(edge))
348  {
349  int e = mapOfEdges.FindIndex(edge);
350  edgeloop->edges.push_back(m_curves[e]);
351  edgeloop->edgeo.push_back(exp.Orientation() == TopAbs_FORWARD
354  }
355  exp.Next();
356  }
357  edgeloops.push_back(edgeloop);
358  }
359 
360  int tote = 0;
361  for (int k = 0; k < edgeloops.size(); k++)
362  {
363  tote += edgeloops[k]->edges.size();
364  }
365 
366  ASSERTL0(tote != 1, "cannot handle periodic curves");
367 
368  CADSurf::OrientateEdges(newSurf, edgeloops);
369  newSurf->SetEdges(edgeloops);
370 
371  // now the loops are orientated, tell the curves how they are
372  for (int k = 0; k < edgeloops.size(); k++)
373  {
374  for (int j = 0; j < edgeloops[k]->edges.size(); j++)
375  {
376  edgeloops[k]->edges[j]->SetAdjSurf(
377  make_pair(newSurf, edgeloops[k]->edgeo[j]));
378  }
379  }
380 
381  m_surfs[i] = newSurf;
382 }
383 
384 Array<OneD, NekDouble> CADSystemOCE::GetBoundingBox()
385 {
386  Array<OneD, NekDouble> bound(6);
387  bound[0] = numeric_limits<double>::max(); // xmin
388  bound[1] = numeric_limits<double>::min(); // xmax
389  bound[2] = numeric_limits<double>::max(); // ymin
390  bound[3] = numeric_limits<double>::min(); // ymax
391  bound[4] = numeric_limits<double>::max(); // zmin
392  bound[5] = numeric_limits<double>::min(); // zmax
393 
394  for (int i = 1; i <= m_curves.size(); i++)
395  {
396  CADCurveSharedPtr c = GetCurve(i);
397  Array<OneD, NekDouble> ends = c->GetMinMax();
398 
399  bound[0] = min(bound[0], min(ends[0], ends[3]));
400  bound[1] = max(bound[1], max(ends[0], ends[3]));
401 
402  bound[2] = min(bound[2], min(ends[1], ends[4]));
403  bound[3] = max(bound[3], max(ends[1], ends[4]));
404 
405  bound[4] = min(bound[4], min(ends[2], ends[5]));
406  bound[5] = max(bound[5], max(ends[2], ends[5]));
407  }
408 
409  return bound;
410 }
411 
412 TopoDS_Shape CADSystemOCE::BuildNACA(string naca)
413 {
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");
418 
419  int n = std::stoi(naca);
420  NekDouble T = (n % 100) / 100.0;
421  n /= 100;
422  NekDouble P = (n % 10) / 10.0;
423  n /= 10;
424  NekDouble M = (n % 10) / 100.0;
425 
426  int np = 25;
427 
428  Array<OneD, NekDouble> xc(np);
429  NekDouble dtheta = M_PI / (np - 1);
430  for (int i = 0; i < np; i++)
431  {
432  xc[i] = (1.0 - cos(i * dtheta)) / 2.0;
433  }
434 
435  Array<OneD, NekDouble> yc(np), dyc(np);
436  for (int i = 0; i < np; i++)
437  {
438  if (xc[i] < P)
439  {
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]);
442  }
443  else
444  {
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]);
448  }
449  }
450 
451  Array<OneD, NekDouble> yt(np);
452  for (int i = 0; i < np; i++)
453  {
454  yt[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]);
458  }
459 
460  Array<OneD, NekDouble> x(2 * np - 1), y(2 * np - 1);
461  int l = 0;
462  for (int i = np - 1; i >= 0; i--, l++)
463  {
464  NekDouble theta = atan(dyc[i]);
465 
466  x[l] = xc[i] - yt[i] * sin(theta);
467  y[l] = yc[i] + yt[i] * cos(theta);
468  }
469  for (int i = 1; i < np; i++)
470  {
471  NekDouble theta = atan(dyc[i]);
472 
473  x[i + np - 1] = xc[i] + yt[i] * sin(theta);
474  y[i + np - 1] = yc[i] - yt[i] * cos(theta);
475  }
476 
477  TColgp_Array1OfPnt pointArray(0, 2 * np - 2);
478 
479  for (int i = 0; i < 2 * np - 1; i++)
480  {
481  pointArray.SetValue(i, gp_Pnt(x[i] * 1000.0, y[i] * 1000.0, 0.0));
482  }
483 
484  GeomAPI_PointsToBSpline spline(pointArray);
485  Handle(Geom_BSplineCurve) curve = spline.Curve();
486 
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();
493 
494  BRepBuilderAPI_MakeWire aeroWireBuilder(aeroEdge, TeEdge);
495  TopoDS_Wire aeroWire = aeroWireBuilder.Wire();
496 
497  gp_Trsf transform;
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);
501  aeroWire.Move(mv);
502 
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();
507 
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();
512 
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();
517 
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();
522 
523  BRepBuilderAPI_MakeWire domWireBuilder(inlEdge, topEdge, outEdge, botEdge);
524  TopoDS_Wire domWire = domWireBuilder.Wire();
525 
526  BRepBuilderAPI_MakeFace face(domWire, true);
527  face.Add(aeroWire);
528 
529  ShapeFix_Face sf(face.Face());
530  sf.FixOrientation();
531 
532  return sf.Face();
533 }
534 
535 TopoDS_Shape CADSystemOCE::BuildGeo(string geo)
536 {
537  ifstream f;
538  f.open(geo.c_str());
539 
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;
548 
549  string fline;
550  string flinetmp;
551 
552  while (!f.eof())
553  {
554  getline(f, fline);
555 
556  boost::erase_all(fline, "\r");
557 
558  vector<string> tmp1, tmp2;
559  boost::split(tmp1, fline, boost::is_any_of("//"));
560  fline = tmp1[0];
561 
562  if (!boost::contains(fline, ";"))
563  {
564  flinetmp += fline;
565  continue;
566  }
567 
568  fline = flinetmp + fline;
569  flinetmp.clear();
570 
571  boost::split(tmp1, fline, boost::is_any_of("="));
572 
573  boost::split(tmp2, tmp1[0], boost::is_any_of("("));
574 
575  string type = tmp2[0];
576  boost::erase_all(tmp2[1], ")");
577  boost::erase_all(tmp2[1], " ");
578  int id = std::stoi(tmp2[1]);
579 
580  boost::erase_all(tmp1[1], " ");
581  boost::erase_all(tmp1[1], "{");
582  boost::erase_all(tmp1[1], "}");
583  boost::erase_all(tmp1[1], ";");
584 
585  string var = tmp1[1];
586 
587  if (boost::iequals(type, "Point"))
588  {
589  points[id] = var;
590  }
591  else if (boost::iequals(type, "Line"))
592  {
593  lines[id] = var;
594  }
595  else if (boost::iequals(type, "Spline"))
596  {
597  splines[id] = var;
598  }
599  else if (boost::iequals(type, "BSpline"))
600  {
601  bsplines[id] = var;
602  }
603  else if (boost::iequals(type, "Circle"))
604  {
605  circles[id] = var;
606  }
607  else if (boost::iequals(type, "Ellipse"))
608  {
609  ellipses[id] = var;
610  }
611  else if (boost::iequals(type, "Line Loop"))
612  {
613  // line loops sometimes have negative entries for gmsh
614  // orientaton purposes
615  // we dont care so remove it
616  boost::erase_all(var, "-");
617  loops[id] = var;
618  }
619  else if (boost::iequals(type, "Plane Surface"))
620  {
621  surfs[id] = var;
622  }
623  else
624  {
625  cout << "not sure " << type << endl;
626  }
627  }
628 
629  map<int, string>::iterator it;
630 
631  // build points
632  map<int, gp_Pnt> cPoints;
633  for (it = points.begin(); it != points.end(); it++)
634  {
635  vector<NekDouble> data;
636  ParseUtils::GenerateVector(it->second, data);
637 
638  cPoints[it->first] =
639  gp_Pnt(data[0] * 1000.0, data[1] * 1000.0, data[2] * 1000.0);
640  }
641 
642  // build edges
643  map<int, TopoDS_Edge> cEdges;
644  for (it = lines.begin(); it != lines.end(); it++)
645  {
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();
650  }
651  for (it = splines.begin(); it != splines.end(); it++)
652  {
653  vector<unsigned int> data;
654  ParseUtils::GenerateVector(it->second, data);
655 
656  TColgp_Array1OfPnt pointArray(0, data.size() - 1);
657 
658  for (int i = 0; i < data.size(); i++)
659  {
660  pointArray.SetValue(i, cPoints[data[i]]);
661  }
662  GeomAPI_PointsToBSpline spline(pointArray);
663  Handle(Geom_BSplineCurve) curve = spline.Curve();
664 
665  BRepBuilderAPI_MakeEdge em(curve);
666  cEdges[it->first] = em.Edge();
667  }
668  for (it = bsplines.begin(); it != bsplines.end(); it++)
669  {
670  vector<unsigned int> data;
671  ParseUtils::GenerateVector(it->second, data);
672 
673  TColgp_Array1OfPnt pointArray(0, data.size() - 1);
674 
675  for (int i = 0; i < data.size(); i++)
676  {
677  pointArray.SetValue(i, cPoints[data[i]]);
678  }
679  Handle(Geom_BezierCurve) curve = new Geom_BezierCurve(pointArray);
680 
681  BRepBuilderAPI_MakeEdge em(curve);
682  cEdges[it->first] = em.Edge();
683  }
684  for (it = circles.begin(); it != circles.end(); it++)
685  {
686  vector<unsigned int> data;
687  ParseUtils::GenerateVector(it->second, data);
688 
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]];
693 
694  NekDouble r1 = start.Distance(centre);
695  NekDouble r2 = end.Distance(centre);
696  ASSERTL0(fabs(r1 - r2) < 1e-7, "Non-matching radii");
697 
698  gp_Circ c;
699  c.SetLocation(centre);
700  c.SetRadius(r1);
701  Handle(Geom_Circle) gc = new Geom_Circle(c);
702 
703  ShapeAnalysis_Curve sac;
704  NekDouble p1, p2;
705  sac.Project(gc, start, 1e-8, start, p1);
706  sac.Project(gc, end, 1e-8, end, p2);
707 
708  // Make sure the arc is always of length less than pi
709  if ((p1 > p2) ^ (fabs(p2 - p1) > M_PI))
710  {
711  swap(p1, p2);
712  }
713 
714  Handle(Geom_TrimmedCurve) tc = new Geom_TrimmedCurve(gc, p1, p2, false);
715 
716  BRepBuilderAPI_MakeEdge em(tc);
717  em.Build();
718  cEdges[it->first] = em.Edge();
719  }
720  for (it = ellipses.begin(); it != ellipses.end(); it++)
721  {
722  vector<unsigned int> data;
723  ParseUtils::GenerateVector(it->second, data);
724 
725  ASSERTL0(data.size() == 4, "Wrong definition of ellipse arc");
726  gp_Pnt start = cPoints[data[0]];
727  gp_Pnt centre = cPoints[data[1]];
728  // data[2] useless??
729  gp_Pnt end = cPoints[data[3]];
730 
731  NekDouble major = start.Distance(centre);
732 
733  gp_Vec v1(centre, start);
734  gp_Vec v2(centre, end);
735 
736  gp_Vec vx(1.0, 0.0, 0.0);
737  NekDouble angle = v1.Angle(vx);
738  // Check for negative rotation
739  if (v1.Y() < 0)
740  {
741  angle *= -1;
742  }
743 
744  v2.Rotate(gp_Ax1(), angle);
745  NekDouble minor = fabs(v2.Y() / sqrt(1.0 - v2.X() * v2.X() / (major * major)));
746 
747  gp_Elips e;
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);
753 
754  ShapeAnalysis_Curve sac;
755  NekDouble p1, p2;
756  sac.Project(ge, start, 1e-8, start, p1);
757  sac.Project(ge, end, 1e-8, end, p2);
758 
759  // Make sure the arc is always of length less than pi
760  if (fabs(p2 - p1) > M_PI)
761  {
762  swap(p1, p2);
763  }
764 
765  Handle(Geom_TrimmedCurve) tc = new Geom_TrimmedCurve(ge, p1, p2, true);
766 
767  BRepBuilderAPI_MakeEdge em(tc);
768  em.Build();
769  cEdges[it->first] = em.Edge();
770  }
771 
772  // build wires
773  map<int, TopoDS_Wire> cWires;
774  for (it = loops.begin(); it != loops.end(); it++)
775  {
776  vector<unsigned int> data;
777  ParseUtils::GenerateVector(it->second, data);
778  BRepBuilderAPI_MakeWire wm;
779  for (int i = 0; i < data.size(); i++)
780  {
781  wm.Add(cEdges[data[i]]);
782  }
783  cWires[it->first] = wm.Wire();
784  }
785 
786  // make surface, at this point assuming its 2D (therefore only 1)
787  // also going to assume that the first loop in the list is the outer domain
788  ASSERTL0(surfs.size() == 1, "more than 1 surf");
789  it = surfs.begin();
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++)
794  {
795  face.Add(cWires[data[i]]);
796  }
797 
798  ASSERTL0(face.Error() == BRepBuilderAPI_FaceDone, "build geo failed");
799 
800  ShapeFix_Face sf(face.Face());
801  sf.FixOrientation();
802 
803  return sf.Face();
804 }
805 }
806 }
CADCurveFactory & GetCADCurveFactory()
Definition: CADSystem.cpp:59
std::shared_ptr< CADSurf > CADSurfSharedPtr
Definition: CADCurve.h:51
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
std::shared_ptr< EdgeLoop > EdgeLoopSharedPtr
Definition: CADSurf.h:65
std::shared_ptr< CADVert > CADVertSharedPtr
Definition: CADCurve.h:49
STL namespace.
struct which descibes a collection of cad edges which are a loop on the cad surface ...
Definition: CADSurf.h:57
CADVertFactory & GetCADVertFactory()
Definition: CADSystem.cpp:53
void filterModShape(TopTools_DataMapOfShapeShape &modShape, TopoDS_Shape &S)
CADSurfFactory & GetCADSurfFactory()
Definition: CADSystem.cpp:65
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
std::shared_ptr< CADCurve > CADCurveSharedPtr
Definition: CADCurve.h:219
double NekDouble
EngineFactory & GetEngineFactory()
Definition: CADSystem.cpp:47
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:199