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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: cad object methods.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
36 
42 
43 using namespace std;
44 
45 namespace Nektar
46 {
47 namespace NekMeshUtils
48 {
49 
50 std::string CADSystemOCE::key = GetEngineFactory().RegisterCreatorFunction(
51  "oce", CADSystemOCE::create, "Uses OCE as cad engine");
52 
53 bool CADSystemOCE::LoadCAD()
54 {
55  if (m_naca.size() == 0)
56  {
57  // not a naca profile behave normally
58  // Takes step file and makes OpenCascade shape
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();
65  if (shape.IsNull())
66  {
67  return false;
68  }
69  }
70  else
71  {
72  shape = BuildNACA(m_name);
73  }
74 
75  TopExp_Explorer explr;
76 
77  // build map of verticies
78  for (explr.Init(shape, TopAbs_VERTEX); explr.More(); explr.Next())
79  {
80  TopoDS_Shape v = explr.Current();
81  if (mapOfVerts.Contains(v))
82  {
83  continue;
84  }
85  int i = mapOfVerts.Add(v);
86  AddVert(i, v);
87  }
88 
89  // For each face of the geometry, get the local edges which bound it. If
90  // they are valid (their type != 7), then add them to an edge map. This
91  // filters out the dummy edges which OCC uses.
92  for (explr.Init(shape, TopAbs_EDGE); explr.More(); explr.Next())
93  {
94  TopoDS_Shape e = explr.Current().Oriented(TopAbs_FORWARD);
95  if (mapOfEdges.Contains(e))
96  {
97  continue;
98  }
99  BRepAdaptor_Curve curve = BRepAdaptor_Curve(TopoDS::Edge(e));
100  if (curve.GetType() != 7)
101  {
102  int i = mapOfEdges.Add(e);
103  AddCurve(i, e);
104  }
105  }
106 
107  for (explr.Init(shape, TopAbs_FACE); explr.More(); explr.Next())
108  {
109  TopoDS_Shape f = explr.Current();
110  ASSERTL0(!mapOfFaces.Contains(f), "duplicated faces");
111  int i = mapOfFaces.Add(f);
112 
113  AddSurf(i, f);
114  }
115 
116  // attempts to identify properties of the vertex on the degen edge
117  for (int i = 1; i <= mapOfFaces.Extent(); i++)
118  {
119  TopoDS_Shape face = mapOfFaces.FindKey(i).Oriented(TopAbs_FORWARD);
120 
121  TopTools_IndexedMapOfShape localEdges;
122  TopExp::MapShapes(face, TopAbs_EDGE, localEdges);
123 
124  for (int j = 1; j <= localEdges.Extent(); j++)
125  {
126  TopoDS_Shape edge = localEdges.FindKey(j);
127  if (BRep_Tool::Degenerated(TopoDS::Edge(edge)))
128  {
129  gp_Pnt2d p1, p2;
130 
131  BRep_Tool::UVPoints(TopoDS::Edge(edge), TopoDS::Face(face), p1,
132  p2);
133 
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);
138  }
139  }
140  }
141 
142  // This checks that all edges are bound by two surfaces, sanity check.
143  if (!m_2d)
144  {
146  for (it = m_curves.begin(); it != m_curves.end(); it++)
147  {
148  ASSERTL0(it->second->GetAdjSurf().size() == 2,
149  "curve is not joined to 2 surfaces");
150  }
151  }
152 
153  return true;
154 }
155 
156 void CADSystemOCE::AddVert(int i, TopoDS_Shape in)
157 {
159 
160  boost::static_pointer_cast<CADVertOCE>(newVert)->Initialise(i, in);
161 
162  m_verts[i] = newVert;
163 }
164 
165 void CADSystemOCE::AddCurve(int i, TopoDS_Shape in)
166 {
168  boost::static_pointer_cast<CADCurveOCE>(newCurve)->Initialise(i, in);
169 
170  TopoDS_Vertex fv = TopExp::FirstVertex(TopoDS::Edge(in));
171  TopoDS_Vertex lv = TopExp::LastVertex(TopoDS::Edge(in));
172 
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);
177 
178  m_curves[i] = newCurve;
179 }
180 
181 void CADSystemOCE::AddSurf(int i, TopoDS_Shape in)
182 {
184  boost::static_pointer_cast<CADSurfOCE>(newSurf)->Initialise(i, in);
185 
186  // do the exploration on forward oriented
187  TopoDS_Shape face = in.Oriented(TopAbs_FORWARD);
188  TopTools_IndexedMapOfShape mapOfWires;
189  TopExp::MapShapes(face, TopAbs_WIRE, mapOfWires);
190  vector<EdgeLoopSharedPtr> edgeloops;
191  // now we acutally analyse the loops for cad building
192  for (int j = 1; j <= mapOfWires.Extent(); j++)
193  {
195 
196  TopoDS_Shape wire = mapOfWires.FindKey(j);
197 
198  BRepTools_WireExplorer exp;
199 
200  exp.Init(TopoDS::Wire(wire));
201 
202  while (exp.More())
203  {
204  TopoDS_Shape edge = exp.Current();
205 
206  if (mapOfEdges.Contains(edge))
207  {
208  int e = mapOfEdges.FindIndex(edge);
209  edgeloop->edges.push_back(m_curves[e]);
210  edgeloop->edgeo.push_back(exp.Orientation() == TopAbs_FORWARD
213  }
214  exp.Next();
215  }
216  edgeloops.push_back(edgeloop);
217  }
218 
219  int tote = 0;
220  for (int k = 0; k < edgeloops.size(); k++)
221  {
222  tote += edgeloops[k]->edges.size();
223  }
224 
225  ASSERTL0(tote != 1, "cannot handle periodic curves");
226 
227  CADSurf::OrientateEdges(newSurf, edgeloops);
228  newSurf->SetEdges(edgeloops);
229 
230  // now the loops are orientated, tell the curves how they are
231  for (int k = 0; k < edgeloops.size(); k++)
232  {
233  for (int j = 0; j < edgeloops[k]->edges.size(); j++)
234  {
235  edgeloops[k]->edges[j]->SetAdjSurf(
236  make_pair(newSurf, edgeloops[k]->edgeo[j]));
237  }
238  }
239 
240  m_surfs[i] = newSurf;
241 }
242 
243 Array<OneD, NekDouble> CADSystemOCE::GetBoundingBox()
244 {
245  Array<OneD, NekDouble> bound(6);
246  bound[0] = numeric_limits<double>::max(); // xmin
247  bound[1] = numeric_limits<double>::min(); // xmax
248  bound[2] = numeric_limits<double>::max(); // ymin
249  bound[3] = numeric_limits<double>::min(); // ymax
250  bound[4] = numeric_limits<double>::max(); // zmin
251  bound[5] = numeric_limits<double>::min(); // zmax
252 
253  for (int i = 1; i <= m_curves.size(); i++)
254  {
255  CADCurveSharedPtr c = GetCurve(i);
256  Array<OneD, NekDouble> ends = c->GetMinMax();
257 
258  bound[0] = min(bound[0], min(ends[0], ends[3]));
259  bound[1] = max(bound[1], max(ends[0], ends[3]));
260 
261  bound[2] = min(bound[2], min(ends[1], ends[4]));
262  bound[3] = max(bound[3], max(ends[1], ends[4]));
263 
264  bound[4] = min(bound[4], min(ends[2], ends[5]));
265  bound[5] = max(bound[5], max(ends[2], ends[5]));
266  }
267 
268  return bound;
269 }
270 
271 TopoDS_Shape CADSystemOCE::BuildNACA(string naca)
272 {
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");
277 
278  int n = boost::lexical_cast<int>(naca);
279  NekDouble T = (n % 100) / 100.0;
280  n /= 100;
281  NekDouble P = (n % 10) / 10.0;
282  n /= 10;
283  NekDouble M = (n % 10) / 100.0;
284 
285  int np = 25;
286 
287  Array<OneD, NekDouble> xc(np);
288  NekDouble dtheta = M_PI / (np - 1);
289  for (int i = 0; i < np; i++)
290  {
291  xc[i] = (1.0 - cos(i * dtheta)) / 2.0;
292  }
293 
294  Array<OneD, NekDouble> yc(np), dyc(np);
295  for (int i = 0; i < np; i++)
296  {
297  if (xc[i] < P)
298  {
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]);
301  }
302  else
303  {
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]);
307  }
308  }
309 
310  Array<OneD, NekDouble> yt(np);
311  for (int i = 0; i < np; i++)
312  {
313  yt[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]);
317  }
318 
319  Array<OneD, NekDouble> x(2 * np - 1), y(2 * np - 1);
320  int l = 0;
321  for (int i = np - 1; i >= 0; i--, l++)
322  {
323  NekDouble theta = atan(dyc[i]);
324 
325  x[l] = xc[i] - yt[i] * sin(theta);
326  y[l] = yc[i] + yt[i] * cos(theta);
327  }
328  for (int i = 1; i < np; i++)
329  {
330  NekDouble theta = atan(dyc[i]);
331 
332  x[i + np - 1] = xc[i] + yt[i] * sin(theta);
333  y[i + np - 1] = yc[i] - yt[i] * cos(theta);
334  }
335 
336  TColgp_Array1OfPnt pointArray(0, 2 * np - 2);
337 
338  for (int i = 0; i < 2 * np - 1; i++)
339  {
340  pointArray.SetValue(i, gp_Pnt(x[i] * 1000.0, y[i] * 1000.0, 0.0));
341  }
342 
343  GeomAPI_PointsToBSpline spline(pointArray);
344  Handle(Geom_BSplineCurve) curve = spline.Curve();
345 
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();
352 
353  BRepBuilderAPI_MakeWire aeroWireBuilder(aeroEdge, TeEdge);
354  TopoDS_Wire aeroWire = aeroWireBuilder.Wire();
355 
356  gp_Trsf transform;
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);
360  aeroWire.Move(mv);
361 
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();
366 
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();
371 
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();
376 
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();
381 
382  BRepBuilderAPI_MakeWire domWireBuilder(inlEdge, topEdge, outEdge, botEdge);
383  TopoDS_Wire domWire = domWireBuilder.Wire();
384 
385  BRepBuilderAPI_MakeFace face(domWire, true);
386  face.Add(aeroWire);
387 
388  ShapeFix_Face sf(face.Face());
389  sf.FixOrientation();
390 
391 
392  return sf.Face();
393 }
394 }
395 }
CADCurveFactory & GetCADCurveFactory()
Definition: CADSystem.cpp:64
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
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.
Definition: NekFactory.hpp:162
STL namespace.
CADVertFactory & GetCADVertFactory()
Definition: CADSystem.cpp:56
CADSurfFactory & GetCADSurfFactory()
Definition: CADSystem.cpp:72
double NekDouble
boost::shared_ptr< CADSurf > CADSurfSharedPtr
Definition: CADSurf.h:172
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
EngineFactory & GetEngineFactory()
Definition: CADSystem.cpp:48
struct which descibes a collection of cad edges which are a loop on the cad surface ...
Definition: CADSystem.h:78
boost::shared_ptr< EdgeLoop > EdgeLoopSharedPtr
Definition: CADSystem.h:86
boost::shared_ptr< CADVert > CADVertSharedPtr
Definition: CADSystem.h:54
boost::shared_ptr< CADCurve > CADCurveSharedPtr
Definition: CADCurve.h:194
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215