Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CADSystem.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 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <boost/algorithm/string.hpp>
37 #include <boost/filesystem.hpp>
38 
44 
45 using namespace std;
46 
47 namespace Nektar
48 {
49 namespace NekMeshUtils
50 {
51 
52 string CADSystem::GetName()
53 {
54  return m_name;
55 }
56 
57 void CADSystem::Report()
58 {
59  cout << endl << "CAD report:" << endl;
60  cout << "\tCAD has: " << m_curves.size() << " curves." << endl;
61  cout << "\tCAD has: " << m_surfs.size() << " surfaces." << endl;
62 }
63 
64 Array<OneD, NekDouble> CADSystem::GetBoundingBox()
65 {
66  Array<OneD, NekDouble> bound(6);
67  bound[0] = numeric_limits<double>::max(); // xmin
68  bound[1] = numeric_limits<double>::min(); // xmax
69  bound[2] = numeric_limits<double>::max(); // ymin
70  bound[3] = numeric_limits<double>::min(); // ymax
71  bound[4] = numeric_limits<double>::max(); // zmin
72  bound[5] = numeric_limits<double>::min(); // zmax
73 
74  for (int i = 1; i <= m_curves.size(); i++)
75  {
76  gp_Pnt start, end;
77  CADCurveSharedPtr c = GetCurve(i);
78  Array<OneD, NekDouble> ends = c->GetMinMax();
79 
80  bound[0] = min(bound[0], min(ends[0], ends[3]));
81  bound[1] = max(bound[1], max(ends[0], ends[3]));
82 
83  bound[2] = min(bound[2], min(ends[1], ends[4]));
84  bound[3] = max(bound[3], max(ends[1], ends[4]));
85 
86  bound[4] = min(bound[4], min(ends[2], ends[5]));
87  bound[5] = max(bound[5], max(ends[2], ends[5]));
88  }
89 
90  return bound;
91 }
92 
93 bool CADSystem::LoadCAD()
94 {
95  if (!boost::filesystem::exists(m_name.c_str()))
96  {
97  return false;
98  }
99 
100  string ext;
101  size_t pos = m_name.find(".");
102  ext = m_name.substr(pos);
103 
104  if (boost::iequals(ext, ".STEP") || boost::iequals(ext, ".STP"))
105  {
106  // Takes step file and makes OpenCascade shape
107  STEPControl_Reader reader;
108  reader = STEPControl_Reader();
109  reader.ReadFile(m_name.c_str());
110  reader.NbRootsForTransfer();
111  reader.TransferRoots();
112  shape = reader.OneShape();
113  if (shape.IsNull())
114  {
115  return false;
116  }
117  }
118  else if (boost::iequals(ext, ".IGES") || boost::iequals(ext, ".IGS"))
119  {
120  // Takes IGES file and makes OpenCascade shape
121  IGESControl_Reader reader;
122  reader = IGESControl_Reader();
123  reader.ReadFile(m_name.c_str());
124  reader.NbRootsForTransfer();
125  reader.TransferRoots();
126  shape = reader.OneShape();
127  if (shape.IsNull())
128  {
129  return false;
130  }
131  }
132  else
133  {
134  return false;
135  }
136 
137  // faces and verts can be extracted straight from shape
138  TopTools_IndexedMapOfShape mapOfVerts, mapOfFaces;
139  TopExp::MapShapes(shape, TopAbs_VERTEX, mapOfVerts);
140  TopExp::MapShapes(shape, TopAbs_FACE, mapOfFaces);
141 
142  // edges need to be built from loops around faces to elimiate degen and
143  // hanging edges
144  TopTools_IndexedMapOfShape mapOfEdges;
145 
146  // build map of verticies
147  for (int i = 1; i <= mapOfVerts.Extent(); i++)
148  {
149  TopoDS_Shape v = mapOfVerts.FindKey(i);
150  AddVert(i, v);
151  }
152 
153  // For each face of the geometry, get the local edges which bound it. If
154  // they are valid (their type != 7), then add them to an edge map. This
155  // filters out the dummy edges which OCC uses.
156  for (int i = 1; i <= mapOfFaces.Extent(); i++)
157  {
158  TopoDS_Shape face = mapOfFaces.FindKey(i);
159 
160  TopTools_IndexedMapOfShape localEdges;
161  TopExp::MapShapes(face, TopAbs_EDGE, localEdges);
162 
163  for (int j = 1; j <= localEdges.Extent(); j++)
164  {
165  TopoDS_Shape edge = localEdges.FindKey(j);
166  BRepAdaptor_Curve curve = BRepAdaptor_Curve(TopoDS::Edge(edge));
167  if (curve.GetType() != 7)
168  {
169  if (!(mapOfEdges.Contains(edge)))
170  {
171  mapOfEdges.Add(edge);
172  }
173  }
174  }
175  }
176 
177  map<int, vector<int> > adjsurfmap; // from id of curve to list of ids of
178  // surfs
179 
180  // Adds edges to our type and map
181  for (int i = 1; i <= mapOfEdges.Extent(); i++)
182  {
183  TopoDS_Shape edge = mapOfEdges.FindKey(i);
184  TopoDS_Vertex fv =
185  TopExp::FirstVertex(TopoDS::Edge(edge), Standard_True);
186  TopoDS_Vertex lv =
187  TopExp::LastVertex(TopoDS::Edge(edge), Standard_True);
188 
189  if (edge.Orientation() == 0)
190  {
191  AddCurve(i, edge, mapOfVerts.FindIndex(fv),
192  mapOfVerts.FindIndex(lv));
193  }
194  else
195  {
196  AddCurve(i, edge, mapOfVerts.FindIndex(lv),
197  mapOfVerts.FindIndex(fv));
198  }
199  }
200 
201  // For each face, examine all the wires (i.e. bounding loops) and
202  // investigates the loop. Using this information, connectivity is determined
203  // and edges are associated with surfaces.
204  for (int i = 1; i <= mapOfFaces.Extent(); i++)
205  {
206  TopoDS_Shape face = mapOfFaces.FindKey(i);
207 
208  TopTools_IndexedMapOfShape mapOfWires;
209  TopExp::MapShapes(face, TopAbs_WIRE, mapOfWires);
210 
211  // this pice of code does an idiot check on the loops to make sure
212  // they dont cross or touch
213  if (mapOfWires.Extent() > 1)
214  {
215  TopoDS_Wire ow = BRepTools::OuterWire(TopoDS::Face(face));
216 
217  vector<TopoDS_Shape> wirefacecuts;
218  vector<gp_Pnt> centersofcutfaces;
219 
220  for (int j = 1; j <= mapOfWires.Extent(); j++)
221  {
222  TopoDS_Shape wire = mapOfWires.FindKey(j);
223 
224  if (wire != ow)
225  {
226  BRepBuilderAPI_MakeFace build(
227  BRep_Tool::Surface(TopoDS::Face(face)), 1e-7);
228  build.Add(TopoDS::Wire(wire));
229  TopoDS_Shape newface = build.Shape();
230  wirefacecuts.push_back(newface);
231  BRepAdaptor_Surface b =
232  BRepAdaptor_Surface(TopoDS::Face(newface));
233  NekDouble u, v;
234  u = (b.LastUParameter() - b.FirstUParameter()) / 2.0;
235  v = (b.LastVParameter() - b.FirstVParameter()) / 2.0;
236  centersofcutfaces.push_back(b.Value(u, v));
237  }
238  }
239  for (int j = 0; j < wirefacecuts.size(); j++)
240  {
241  for (int k = 0; k < wirefacecuts.size(); k++)
242  {
243  if (j == k)
244  continue;
245 
246  /// TODO fix this test
247  BRepClass_FaceClassifier fc(TopoDS::Face(wirefacecuts[j]),
248  centersofcutfaces[k], 1e-7);
249  // ASSERTL0(fc.State() == 1, "Internal face loops make this
250  // cad impossible to mesh");
251  }
252  }
253  }
254 
255  vector<EdgeLoop> edgeloops;
256 
257  // now we acutally analyse the loops for cad building
258  for (int j = 1; j <= mapOfWires.Extent(); j++)
259  {
260  EdgeLoop edgeloop;
261 
262  TopoDS_Shape wire = mapOfWires.FindKey(j);
263 
264  ShapeAnalysis_Wire wiretest(TopoDS::Wire(wire), TopoDS::Face(face),
265  1E-6);
266 
267  // calculate the center of the wire
268  GProp_GProps massProps;
269  BRepGProp::SurfaceProperties(wire, massProps);
270  gp_Pnt gPt = massProps.CentreOfMass();
271  // alternative locuv methods
272  ShapeAnalysis_Surface sas(BRep_Tool::Surface(TopoDS::Face(face)));
273  sas.SetDomain(
274  BRepAdaptor_Surface(TopoDS::Face(face)).FirstUParameter(),
275  BRepAdaptor_Surface(TopoDS::Face(face)).LastUParameter(),
276  BRepAdaptor_Surface(TopoDS::Face(face)).FirstVParameter(),
277  BRepAdaptor_Surface(TopoDS::Face(face)).LastVParameter());
278  gp_Pnt2d p2 = sas.ValueOfUV(gPt, 1e-7);
279  Array<OneD, NekDouble> cen(2);
280  cen[0] = p2.X();
281  cen[1] = p2.Y();
282  edgeloop.center = cen;
283 
284  BRepTools_WireExplorer exp;
285 
286  exp.Init(TopoDS::Wire(wire));
287 
288  while (exp.More())
289  {
290  TopoDS_Shape edge = exp.Current();
291 
292  if (mapOfEdges.Contains(edge))
293  {
294  int e = mapOfEdges.FindIndex(edge);
295  edgeloop.edges.push_back(m_curves[e]);
296  edgeloop.edgeo.push_back(exp.Orientation());
297  adjsurfmap[e].push_back(i);
298  }
299 
300  exp.Next();
301  }
302 
303  edgeloops.push_back(edgeloop);
304  }
305 
306  AddSurf(i, face, edgeloops);
307  }
308 
309  // attempts to identify properties of the vertex on the degen edge
310  for (int i = 1; i <= mapOfFaces.Extent(); i++)
311  {
312  TopoDS_Shape face = mapOfFaces.FindKey(i);
313 
314  TopTools_IndexedMapOfShape localEdges;
315  TopExp::MapShapes(face, TopAbs_EDGE, localEdges);
316 
317  for (int j = 1; j <= localEdges.Extent(); j++)
318  {
319  TopoDS_Shape edge = localEdges.FindKey(j);
320  if (BRep_Tool::Degenerated(TopoDS::Edge(edge)))
321  {
322  cout << "degen edge on face " << i << endl;
323  gp_Pnt2d p1, p2;
324 
325  BRep_Tool::UVPoints(TopoDS::Edge(edge), TopoDS::Face(face), p1,
326  p2);
327 
328  m_verts[mapOfVerts.FindIndex(TopExp::FirstVertex(
329  TopoDS::Edge(edge), Standard_True))]
330  ->SetDegen(i, m_surfs[i], (p1.X() + p2.X()) / 2.0,
331  (p1.Y() + p2.Y()) / 2.0);
332  }
333  }
334  }
335 
336  // This checks that all edges are bound by two surfaces, sanity check.
337  for (map<int, vector<int> >::iterator it = adjsurfmap.begin();
338  it != adjsurfmap.end(); it++)
339  {
340  ASSERTL0(it->second.size() == 2, "no three curve surfaces");
341  vector<CADSurfSharedPtr> sfs;
342  for (int i = 0; i < it->second.size(); i++)
343  {
344  sfs.push_back(m_surfs[it->second[i]]);
345  }
346  m_curves[it->first]->SetAdjSurf(sfs);
347  }
348  return true;
349 }
350 
351 void CADSystem::AddVert(int i, TopoDS_Shape in)
352 {
354 
355  m_verts[i] = newVert;
356 }
357 
358 void CADSystem::AddCurve(int i, TopoDS_Shape in, int fv, int lv)
359 {
360  CADCurveSharedPtr newCurve =
362 
363  vector<CADVertSharedPtr> vs;
364  vs.push_back(m_verts[fv]);
365  vs.push_back(m_verts[lv]);
366  m_curves[i] = newCurve;
367  m_curves[i]->SetVert(vs);
368 }
369 
370 void CADSystem::AddSurf(int i, TopoDS_Shape in, vector<EdgeLoop> ein)
371 {
372  CADSurfSharedPtr newSurf =
374  m_surfs[i] = newSurf;
375 
376  if (in.Orientation() == 0)
377  {
378  m_surfs[i]->SetReverseNomral();
379  }
380 
381  int tote = 0;
382  for (int i = 0; i < ein.size(); i++)
383  {
384  tote += ein[i].edges.size();
385  }
386 
387  ASSERTL0(tote != 1, "cannot handle periodic curves");
388 
389  if (tote == 2)
390  {
391  m_surfs[i]->SetTwoC();
392  }
393 }
394 
395 bool CADSystem::InsideShape(Array<OneD, NekDouble> loc)
396 {
397  gp_Pnt p(loc[0] * 1000.0, loc[1] * 1000.0, loc[2] * 1000.0);
398 
399  BRepClass3d_SolidClassifier test(shape, p, 1E-7);
400  if (test.State() == TopAbs_IN || test.State() == TopAbs_ON)
401  {
402  return true;
403  }
404  else
405  {
406  return false;
407  }
408 }
409 }
410 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Array< OneD, NekDouble > center
Definition: CADSurf.h:64
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
STL namespace.
struct which descibes a collection of cad edges which for a loop on the cad surface ...
Definition: CADSurf.h:60
std::vector< CADCurveSharedPtr > edges
Definition: CADSurf.h:62
double NekDouble
boost::shared_ptr< CADSurf > CADSurfSharedPtr
Definition: CADSurf.h:217
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< CADVert > CADVertSharedPtr
Definition: CADSystem.h:51
std::vector< int > edgeo
Definition: CADSurf.h:63
boost::shared_ptr< CADCurve > CADCurveSharedPtr
Definition: CADCurve.h:168