Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CADSurf.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 
36 #include <sstream>
37 
39 
40 using namespace std;
41 
42 namespace Nektar
43 {
44 namespace NekMeshUtils
45 {
46 
47 CADSurf::CADSurf(int i, TopoDS_Shape in, vector<EdgeLoop> ein) : m_edges(ein)
48 {
49  // this bit of code changes the units of the cad from mm opencascade
50  // defualt to m
51  gp_Trsf transform;
52  gp_Pnt ori(0.0, 0.0, 0.0);
53  transform.SetScale(ori, 1.0 / 1000.0);
54  TopLoc_Location mv(transform);
55  m_s = BRep_Tool::Surface(TopoDS::Face(in));
56  in.Move(mv);
57  m_occSurface = BRepAdaptor_Surface(TopoDS::Face(in));
58  m_correctNormal = true;
59  m_hasTwoCurves = false;
60  m_id = i;
61  m_type = surf;
62 }
63 
65 {
66  // has to transfer back to mm
67  gp_Pnt loc(p[0] * 1000.0, p[1] * 1000.0, p[2] * 1000.0);
68 
69  GeomAPI_ProjectPointOnSurf projection(
70  loc, m_s, m_occSurface.FirstUParameter(), m_occSurface.LastUParameter(),
71  m_occSurface.FirstVParameter(), m_occSurface.LastVParameter(),
72  Extrema_ExtAlgo_Tree);
73 
75  if (projection.NbPoints() == 0)
76  {
77  // alternative locuv methods
78  ShapeAnalysis_Surface sas(m_s);
79  sas.SetDomain(
80  m_occSurface.FirstUParameter(), m_occSurface.LastUParameter(),
81  m_occSurface.FirstVParameter(), m_occSurface.LastVParameter());
82 
83  gp_Pnt2d p2 = sas.ValueOfUV(loc, 1e-7);
84  uvr[0] = p2.X();
85  uvr[1] = p2.Y();
86 
87  gp_Pnt p3 = sas.Value(p2);
88  ASSERTL0(p3.Distance(loc) < 1e-3, "large locuv distance sas");
89  }
90  else
91  {
92  Quantity_Parameter ui;
93  Quantity_Parameter vi;
94 
95  projection.Parameters(1, ui, vi);
96 
97  uvr[0] = ui;
98  uvr[1] = vi;
99 
100  if (projection.Distance(1) > 1.0)
101  {
102  stringstream ss;
103  cerr << "large locuv distance " << projection.Distance(1) / 1000.0
104  << endl;
105  }
106  }
107 
108  // if the uv returned is slightly off the surface
109  //(which ShapeAnalysis_Surface can do sometimes)
110  if (uvr[0] < m_occSurface.FirstUParameter() ||
111  uvr[0] > m_occSurface.LastUParameter() ||
112  uvr[1] < m_occSurface.FirstVParameter() ||
113  uvr[1] > m_occSurface.LastVParameter())
114  {
115  if (uvr[0] < m_occSurface.FirstUParameter() &&
116  fabs(m_occSurface.FirstUParameter() - uvr[0]) < 1E-6)
117  {
118  uvr[0] = m_occSurface.FirstUParameter();
119  }
120  else if (uvr[0] > m_occSurface.LastUParameter() &&
121  fabs(m_occSurface.LastUParameter() - uvr[0]) < 1E-6)
122  {
123  uvr[0] = m_occSurface.LastUParameter();
124  }
125  else if (uvr[1] < m_occSurface.FirstVParameter() &&
126  fabs(m_occSurface.FirstVParameter() - uvr[1]) < 1E-6)
127  {
128  uvr[1] = m_occSurface.FirstVParameter();
129  }
130  else if (uvr[1] > m_occSurface.LastVParameter() &&
131  fabs(m_occSurface.LastVParameter() - uvr[1]) < 1E-6)
132  {
133  uvr[1] = m_occSurface.LastVParameter();
134  }
135  else
136  {
137  ASSERTL0(false, "Cannot correct locuv");
138  }
139  }
140 
141  return uvr;
142 }
143 
145 {
146  Test(uv);
147 
148  Array<OneD, NekDouble> n = N(uv);
149 
150  // a zero normal occurs at a signularity, CurvaturePoint
151  // cannot be sampled here
152  if (n[0] == 0 && n[1] == 0 && n[2] == 0)
153  {
154  return 0.0;
155  }
156 
157  Array<OneD, NekDouble> r = D2(uv);
158 
159  // metric and curvature tensors
160  NekDouble E = r[3] * r[3] + r[4] * r[4] + r[5] * r[5];
161  NekDouble F = r[3] * r[6] + r[4] * r[7] + r[5] * r[8];
162  NekDouble G = r[6] * r[6] + r[7] * r[7] + r[8] * r[8];
163  NekDouble e = n[0] * r[9] + n[1] * r[10] + n[2] * r[11];
164  NekDouble f = n[0] * r[15] + n[1] * r[16] + n[2] * r[17];
165  NekDouble g = n[0] * r[12] + n[1] * r[13] + n[2] * r[14];
166 
167  // if det is zero cannot invert matrix, R=0 so must skip
168  if (E * G - F * F < 1E-30)
169  {
170  return 0.0;
171  }
172 
173  NekDouble K, H;
174 
175  K = (e * g - f * f) / (E * G - F * F);
176  H = 0.5 * (e * G - 2 * f * F + g * E) / (E * G - F * F);
177 
178  NekDouble kv[2];
179  kv[0] = abs(H + sqrt(H * H - K));
180  kv[1] = abs(H - sqrt(H * H - K));
181 
182  return kv[0] > kv[1] ? kv[0] : kv[1];
183 }
184 
186 {
187  gp_Pnt loc(p[0] * 1000.0, p[1] * 1000.0, p[2] * 1000.0);
188 
189  // alternative locuv methods
190  ShapeAnalysis_Surface sas(m_s);
191  sas.SetDomain(m_occSurface.FirstUParameter(), m_occSurface.LastUParameter(),
192  m_occSurface.FirstVParameter(),
193  m_occSurface.LastVParameter());
194 
195  gp_Pnt2d p2 = sas.ValueOfUV(loc, 1e-7);
196 
197  gp_Pnt p3 = sas.Value(p2);
198 
199  return p3.Distance(loc);
200 }
201 
203 {
204  gp_Pnt loc(tp[0] * 1000.0, tp[1] * 1000.0, tp[2] * 1000.0);
205 
206  // alternative locuv methods
207  ShapeAnalysis_Surface sas(m_s);
208  sas.SetDomain(m_occSurface.FirstUParameter(), m_occSurface.LastUParameter(),
209  m_occSurface.FirstVParameter(),
210  m_occSurface.LastVParameter());
211 
212  gp_Pnt2d p2 = sas.ValueOfUV(loc, 1e-7);
213 
214  gp_Pnt p3 = sas.Value(p2);
215 
216  tp[0] = p3.X() / 1000.0;
217  tp[1] = p3.Y() / 1000.0;
218  tp[2] = p3.Z() / 1000.0;
219 
220  uv[0] = p2.X();
221  uv[1] = p2.Y();
222 }
223 
225 {
226  if (m_occSurface.GetType() == GeomAbs_Plane)
227  {
228  return true;
229  }
230  else
231  {
232  return false;
233  }
234 }
235 
237 {
238  Test(uv);
239 
240  Array<OneD, NekDouble> location(3);
241  gp_Pnt loc;
242  loc = m_occSurface.Value(uv[0], uv[1]);
243  location[0] = loc.X();
244  location[1] = loc.Y();
245  location[2] = loc.Z();
246  return location;
247 }
248 
250 {
251  Test(uv);
252 
253  Array<OneD, NekDouble> normal(3);
254  gp_Pnt Loc;
255  gp_Vec D1U, D1V;
256  m_occSurface.D1(uv[0], uv[1], Loc, D1U, D1V);
257  gp_Vec n = D1U.Crossed(D1V);
258 
259  if (!m_correctNormal)
260  {
261  n.Reverse();
262  }
263 
264  if (n.X() == 0 && n.Y() == 0 && n.Z() == 0)
265  {
266  // Return bad normal
267  normal[0] = 0.0;
268  normal[1] = 0.0;
269  normal[2] = 0.0;
270  }
271  else
272  {
273  n.Normalize();
274  normal[0] = n.X();
275  normal[1] = n.Y();
276  normal[2] = n.Z();
277  }
278 
279  return normal;
280 }
281 
283 {
284  Test(uv);
285 
287  gp_Pnt Loc;
288  gp_Vec D1U, D1V;
289  m_occSurface.D1(uv[0], uv[1], Loc, D1U, D1V);
290 
291  r[0] = Loc.X(); // x
292  r[1] = Loc.Y(); // y
293  r[2] = Loc.Z(); // z
294  r[3] = D1U.X(); // dx/du
295  r[4] = D1U.Y(); // dy/du
296  r[5] = D1U.Z(); // dz/du
297  r[6] = D1V.X(); // dx/dv
298  r[7] = D1V.Y(); // dy/dv
299  r[8] = D1V.Z(); // dz/dv
300 
301  return r;
302 }
303 
305 {
306  Test(uv);
307 
309  gp_Pnt Loc;
310  gp_Vec D1U, D1V, D2U, D2V, D2UV;
311  m_occSurface.D2(uv[0], uv[1], Loc, D1U, D1V, D2U, D2V, D2UV);
312 
313  r[0] = Loc.X(); // x
314  r[1] = Loc.Y(); // y
315  r[2] = Loc.Z(); // z
316  r[3] = D1U.X(); // dx/dx
317  r[4] = D1U.Y(); // dy/dy
318  r[5] = D1U.Z(); // dz/dz
319  r[6] = D1V.X(); // dx/dx
320  r[7] = D1V.Y(); // dy/dy
321  r[8] = D1V.Z(); // dz/dz
322  r[9] = D2U.X(); // d2x/du2
323  r[10] = D2U.Y(); // d2y/du2
324  r[11] = D2U.Z(); // d2z/du2
325  r[12] = D2V.X(); // d2x/dv2
326  r[13] = D2V.Y(); // d2y/dv2
327  r[14] = D2V.Z(); // d2z/dv2
328  r[15] = D2UV.X(); // d2x/dudv
329  r[16] = D2UV.Y(); // d2y/dudv
330  r[17] = D2UV.Z(); // d2z/dudv
331 
332  return r;
333 }
334 
336 {
337  stringstream error;
338 
339  error << "Point not within parameter plane: ";
340 
341  bool passed = true;
342 
343  if (uv[0] < m_occSurface.FirstUParameter())
344  {
345  if (fabs(uv[0] - m_occSurface.FirstUParameter()) > 1E-8)
346  {
347  error << "U(" << uv[0] << ") is less than Umin("
348  << m_occSurface.FirstUParameter() << ")";
349  passed = false;
350  }
351  }
352  else if (uv[0] > m_occSurface.LastUParameter())
353  {
354  if (fabs(uv[0] - m_occSurface.LastUParameter()) > 1E-8)
355  {
356  error << "U(" << uv[0] << ") is greater than Umax("
357  << m_occSurface.LastUParameter() << ")";
358  passed = false;
359  }
360  }
361  else if (uv[1] < m_occSurface.FirstVParameter())
362  {
363  if (fabs(uv[1] - m_occSurface.FirstVParameter()) > 1E-8)
364  {
365  error << "V(" << uv[1] << ") is less than Vmin("
366  << m_occSurface.FirstVParameter() << ")";
367  passed = false;
368  }
369  }
370  else if (uv[1] > m_occSurface.LastVParameter())
371  {
372  if (fabs(uv[1] - m_occSurface.LastVParameter()) > 1E-8)
373  {
374  error << "V(" << uv[1] << ") is greater than Vmax("
375  << m_occSurface.LastVParameter() << ")";
376  passed = false;
377  }
378  }
379 
380  error << " On Surface: " << GetId();
381 
382  ASSERTL0(passed, error.str());
383 }
384 }
385 }
Array< OneD, NekDouble > D1(Array< OneD, NekDouble > uv)
Get the set of first derivatives at parametric point u,v.
Definition: CADSurf.cpp:282
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Array< OneD, NekDouble > P(Array< OneD, NekDouble > uv)
Get the x,y,z at parametric point u,v.
Definition: CADSurf.cpp:236
bool m_hasTwoCurves
flag to alert the mesh generation to a potential problem is both curves have only two points in the m...
Definition: CADSurf.h:208
Array< OneD, NekDouble > locuv(Array< OneD, NekDouble > p)
Performs a reverse look up to find u,v and x,y,z.
Definition: CADSurf.cpp:64
STL namespace.
BRepAdaptor_Surface m_occSurface
OpenCascade object for surface.
Definition: CADSurf.h:210
Array< OneD, NekDouble > N(Array< OneD, NekDouble > uv)
Get the normal vector at parametric point u,v.
Definition: CADSurf.cpp:249
NekDouble DistanceTo(Array< OneD, NekDouble > p)
does unconstrained locuv to project point from anywhere
Definition: CADSurf.cpp:185
double NekDouble
int m_id
ID of the vert.
Definition: CADObj.h:91
int GetId()
Return ID of the vertex.
Definition: CADObj.h:79
bool IsPlane()
returns true if the surface is flat (2D)
Definition: CADSurf.cpp:224
void ProjectTo(Array< OneD, NekDouble > &tp, Array< OneD, NekDouble > &uv)
Definition: CADSurf.cpp:202
void Test(Array< OneD, NekDouble > uv)
Function which tests the the value of uv used is within the surface.
Definition: CADSurf.cpp:335
cadType m_type
type of the cad object
Definition: CADObj.h:93
Array< OneD, NekDouble > D2(Array< OneD, NekDouble > uv)
Get the set of second derivatives at parametric point u,v.
Definition: CADSurf.cpp:304
NekDouble Curvature(Array< OneD, NekDouble > uv)
returns curvature at point uv
Definition: CADSurf.cpp:144