Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CADSurfOCE.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: CADSurf.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 surface methods.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
37 
38 using namespace std;
39 
40 namespace Nektar
41 {
42 namespace NekMeshUtils
43 {
44 
45 std::string CADSurfOCE::key = GetCADSurfFactory().RegisterCreatorFunction(
46  "oce", CADSurfOCE::create, "CADSurfOCE");
47 
48 void CADSurfOCE::Initialise(int i, TopoDS_Shape in)
49 {
50  // this bit of code changes the units of the cad from mm opencascade
51  // defualt to m
52 
53  m_s = BRep_Tool::Surface(TopoDS::Face(in));
54 
55  if (in.Orientation() == 1)
56  {
57  m_orientation = CADOrientation::eBackwards;
58  }
59 
60  gp_Trsf transform;
61  gp_Pnt ori(0.0, 0.0, 0.0);
62  transform.SetScale(ori, 1.0 / 1000.0);
63  TopLoc_Location mv(transform);
64 
65  in.Move(mv);
66  m_occSurface = BRepAdaptor_Surface(TopoDS::Face(in));
67  m_id = i;
68 
69  m_bounds = Array<OneD, NekDouble>(4);
70  BRepTools::UVBounds(TopoDS::Face(in), m_bounds[0], m_bounds[1], m_bounds[2],
71  m_bounds[3]);
72  m_sas = new ShapeAnalysis_Surface(m_s);
73  m_sas->SetDomain(m_bounds[0], m_bounds[1], m_bounds[2], m_bounds[3]);
74 }
75 
76 Array<OneD, NekDouble> CADSurfOCE::GetBounds()
77 {
78  return m_bounds;
79 }
80 
82 {
83  // has to transfer back to mm
84  gp_Pnt loc(p[0] * 1000.0, p[1] * 1000.0, p[2] * 1000.0);
85 
87 
88  gp_Pnt2d p2 = m_sas->ValueOfUV(loc, 1e-3);
89  uvr[0] = p2.X();
90  uvr[1] = p2.Y();
91 
92  gp_Pnt p3 = m_sas->Value(p2);
93  WARNINGL2(p3.Distance(loc) < 1e-3, "large locuv distance " +
94  boost::lexical_cast<string>(p3.Distance(loc)/1000.0) + " " +
95  boost::lexical_cast<string>(m_id));
96 
97  // if the uv returned is slightly off the surface
98  //(which ShapeAnalysis_Surface can do sometimes)
99  if (uvr[0] < m_bounds[0] || uvr[0] > m_bounds[1] || uvr[1] < m_bounds[2] ||
100  uvr[1] > m_bounds[3])
101  {
102  if (uvr[0] < m_bounds[0])
103  {
104  uvr[0] = m_bounds[0];
105  }
106  else if (uvr[0] > m_bounds[1])
107  {
108  uvr[0] = m_bounds[1];
109  }
110  else if (uvr[1] < m_bounds[2])
111  {
112  uvr[1] = m_bounds[2];
113  }
114  else if (uvr[1] > m_bounds[3])
115  {
116  uvr[1] = m_bounds[3];
117  }
118  else
119  {
120  ASSERTL0(false, "Cannot correct locuv");
121  }
122  }
123 
124  return uvr;
125 }
126 
127 NekDouble CADSurfOCE::Curvature(Array<OneD, NekDouble> uv)
128 {
129 #if defined(NEKTAR_DEBUG)
130  Test(uv);
131 #endif
132 
133  Array<OneD, NekDouble> n = N(uv);
134 
135  // a zero normal occurs at a signularity, CurvaturePoint
136  // cannot be sampled here
137  if (n[0] == 0 && n[1] == 0 && n[2] == 0)
138  {
139  return 0.0;
140  }
141 
142  Array<OneD, NekDouble> r = D2(uv);
143 
144  // metric and curvature tensors
145  NekDouble E = r[3] * r[3] + r[4] * r[4] + r[5] * r[5];
146  NekDouble F = r[3] * r[6] + r[4] * r[7] + r[5] * r[8];
147  NekDouble G = r[6] * r[6] + r[7] * r[7] + r[8] * r[8];
148  NekDouble e = n[0] * r[9] + n[1] * r[10] + n[2] * r[11];
149  NekDouble f = n[0] * r[15] + n[1] * r[16] + n[2] * r[17];
150  NekDouble g = n[0] * r[12] + n[1] * r[13] + n[2] * r[14];
151 
152  // if det is zero cannot invert matrix, R=0 so must skip
153  if (E * G - F * F < 1E-30)
154  {
155  return 0.0;
156  }
157 
158  NekDouble K, H;
159 
160  K = (e * g - f * f) / (E * G - F * F);
161  H = 0.5 * (e * G - 2 * f * F + g * E) / (E * G - F * F);
162 
163  NekDouble kv[2];
164  kv[0] = abs(H + sqrt(H * H - K));
165  kv[1] = abs(H - sqrt(H * H - K));
166 
167  return kv[0] > kv[1] ? kv[0] : kv[1];
168 }
169 
170 NekDouble CADSurfOCE::DistanceTo(Array<OneD, NekDouble> p)
171 {
172  gp_Pnt loc(p[0] * 1000.0, p[1] * 1000.0, p[2] * 1000.0);
173 
174  // alternative locuv methods
175  ShapeAnalysis_Surface sas(m_s);
176  sas.SetDomain(m_bounds[0], m_bounds[1], m_bounds[2], m_bounds[3]);
177 
178  gp_Pnt2d p2 = sas.ValueOfUV(loc, 1e-7);
179 
180  gp_Pnt p3 = sas.Value(p2);
181 
182  return p3.Distance(loc);
183 }
184 
185 void CADSurfOCE::ProjectTo(Array<OneD, NekDouble> &tp,
187 {
188  gp_Pnt loc(tp[0] * 1000.0, tp[1] * 1000.0, tp[2] * 1000.0);
189 
190  // alternative locuv methods
191  ShapeAnalysis_Surface sas(m_s);
192  sas.SetDomain(m_bounds[0], m_bounds[1], m_bounds[2], m_bounds[3]);
193 
194  gp_Pnt2d p2 = sas.ValueOfUV(loc, 1e-7);
195 
196  gp_Pnt p3 = sas.Value(p2);
197 
198  tp[0] = p3.X() / 1000.0;
199  tp[1] = p3.Y() / 1000.0;
200  tp[2] = p3.Z() / 1000.0;
201 
202  uv[0] = p2.X();
203  uv[1] = p2.Y();
204 }
205 
207 {
208 #if defined(NEKTAR_DEBUG)
209  Test(uv);
210 #endif
211 
212  Array<OneD, NekDouble> location(3);
213  gp_Pnt loc;
214  loc = m_occSurface.Value(uv[0], uv[1]);
215  location[0] = loc.X();
216  location[1] = loc.Y();
217  location[2] = loc.Z();
218  return location;
219 }
220 
222 {
223 #if defined(NEKTAR_DEBUG)
224  Test(uv);
225 #endif
226 
227  BRepLProp_SLProps slp(m_occSurface, 2, 1e-6);
228  slp.SetParameters(uv[0], uv[1]);
229 
230  if (!slp.IsNormalDefined())
231  {
232  return Array<OneD, NekDouble>(3, 0.0);
233  }
234 
235  gp_Dir d = slp.Normal();
236 
237  Array<OneD, NekDouble> normal(3);
238 
239  if (m_orientation == CADOrientation::eBackwards)
240  {
241  d.Reverse();
242  }
243 
244  normal[0] = d.X();
245  normal[1] = d.Y();
246  normal[2] = d.Z();
247 
248  return normal;
249 }
250 
252 {
253 #if defined(NEKTAR_DEBUG)
254  Test(uv);
255 #endif
256 
258  gp_Pnt Loc;
259  gp_Vec D1U, D1V;
260  m_occSurface.D1(uv[0], uv[1], Loc, D1U, D1V);
261 
262  r[0] = Loc.X(); // x
263  r[1] = Loc.Y(); // y
264  r[2] = Loc.Z(); // z
265  r[3] = D1U.X(); // dx/du
266  r[4] = D1U.Y(); // dy/du
267  r[5] = D1U.Z(); // dz/du
268  r[6] = D1V.X(); // dx/dv
269  r[7] = D1V.Y(); // dy/dv
270  r[8] = D1V.Z(); // dz/dv
271 
272  return r;
273 }
274 
276 {
277 #if defined(NEKTAR_DEBUG)
278  Test(uv);
279 #endif
280 
282  gp_Pnt Loc;
283  gp_Vec D1U, D1V, D2U, D2V, D2UV;
284  m_occSurface.D2(uv[0], uv[1], Loc, D1U, D1V, D2U, D2V, D2UV);
285 
286  r[0] = Loc.X(); // x
287  r[1] = Loc.Y(); // y
288  r[2] = Loc.Z(); // z
289  r[3] = D1U.X(); // dx/dx
290  r[4] = D1U.Y(); // dy/dy
291  r[5] = D1U.Z(); // dz/dz
292  r[6] = D1V.X(); // dx/dx
293  r[7] = D1V.Y(); // dy/dy
294  r[8] = D1V.Z(); // dz/dz
295  r[9] = D2U.X(); // d2x/du2
296  r[10] = D2U.Y(); // d2y/du2
297  r[11] = D2U.Z(); // d2z/du2
298  r[12] = D2V.X(); // d2x/dv2
299  r[13] = D2V.Y(); // d2y/dv2
300  r[14] = D2V.Z(); // d2z/dv2
301  r[15] = D2UV.X(); // d2x/dudv
302  r[16] = D2UV.Y(); // d2y/dudv
303  r[17] = D2UV.Z(); // d2z/dudv
304 
305  return r;
306 }
307 
308 void CADSurfOCE::Test(Array<OneD, NekDouble> uv)
309 {
310  stringstream error;
311 
312  error << "Point not within parameter plane: ";
313 
314  bool passed = true;
315 
316  if (uv[0] < m_bounds[0])
317  {
318  if (fabs(uv[0] - m_bounds[0]) > 1E-6)
319  {
320  error << "U(" << uv[0] << ") is less than Umin(" << m_bounds[0]
321  << ")";
322  passed = false;
323  }
324  }
325  else if (uv[0] > m_bounds[1])
326  {
327  if (fabs(uv[0] - m_bounds[1]) > 1E-6)
328  {
329  error << "U(" << uv[0] << ") is greater than Umax(" << m_bounds[1]
330  << ")";
331  passed = false;
332  }
333  }
334  else if (uv[1] < m_bounds[2])
335  {
336  if (fabs(uv[1] - m_bounds[2]) > 1E-6)
337  {
338  error << "V(" << uv[1] << ") is less than Vmin(" << m_bounds[2]
339  << ")";
340  passed = false;
341  }
342  }
343  else if (uv[1] > m_bounds[3])
344  {
345  if (fabs(uv[1] - m_bounds[3]) > 1E-6)
346  {
347  error << "V(" << uv[1] << ") is greater than Vmax(" << m_bounds[3]
348  << ")";
349  passed = false;
350  }
351  }
352 
353  error << " On Surface: " << GetId();
354  ASSERTL1(passed, "Warning: " + error.str());
355 }
356 }
357 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
STL namespace.
CADSurfFactory & GetCADSurfFactory()
Definition: CADSystem.cpp:72
double NekDouble
#define WARNINGL2(condition, msg)
Definition: ErrorUtil.hpp:251
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215