Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CADCurveOCE.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: CADCurve.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 curve methods.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
37 
38 using namespace std;
39 
40 namespace Nektar
41 {
42 namespace NekMeshUtils
43 {
44 
45 std::string CADCurveOCE::key = GetCADCurveFactory().RegisterCreatorFunction(
46  "oce", CADCurveOCE::create, "CADCurveOCE");
47 
48 NekDouble CADCurveOCE::tAtArcLength(NekDouble s)
49 {
50  NekDouble dt =
51  (m_occCurve.LastParameter() - m_occCurve.FirstParameter()) / (1000);
52  NekDouble t = m_occCurve.FirstParameter();
53 
54  NekDouble len = 0.0;
55 
56  while (len <= s)
57  {
58  gp_Pnt P1, P2;
59  gp_Vec drdt1, drdt2;
60 
61  m_occCurve.D1(t, P1, drdt1);
62  t += dt;
63  m_occCurve.D1(t, P2, drdt2);
64 
65  len += (drdt1.Magnitude() + drdt2.Magnitude()) / 2.0 * dt;
66  }
67 
68  return t - dt;
69 }
70 
71 NekDouble CADCurveOCE::Length(NekDouble ti, NekDouble tf)
72 {
73  Array<OneD, NekDouble> b = GetBounds();
74  Handle(Geom_Curve) NewCurve = new Geom_TrimmedCurve(m_c, ti, tf);
75  TopoDS_Edge NewEdge = BRepBuilderAPI_MakeEdge(NewCurve);
76  GProp_GProps System;
77  BRepGProp::LinearProperties(NewEdge, System);
78  return System.Mass() / 1000.0;
79 }
80 
81 NekDouble CADCurveOCE::loct(Array<OneD, NekDouble> xyz)
82 {
83  NekDouble t = 0.0;
84  Array<OneD, NekDouble> b = GetBounds();
85 
86  gp_Pnt loc(xyz[0]*1000.0, xyz[1]*1000.0, xyz[2]*1000.0);
87 
88  ShapeAnalysis_Curve sac;
89  gp_Pnt p;
90  sac.Project(m_c,loc,1e-8 ,p,t);
91 
92  if(p.Distance(loc) > 1e-5)
93  {
94  cerr << "large loct distance" << endl;
95  }
96  return t;
97 }
98 
100 {
101  Array<OneD, NekDouble> location(3);
102  gp_Pnt loc = m_occCurve.Value(t);
103 
104  location[0] = loc.X();
105  location[1] = loc.Y();
106  location[2] = loc.Z();
107 
108  return location;
109 }
110 
112 {
113  Array<OneD, NekDouble> out(9);
114  gp_Pnt loc;
115  gp_Vec d1, d2;
116  m_occCurve.D2(t, loc, d1, d2);
117 
118  out[0] = loc.X();
119  out[1] = loc.Y();
120  out[2] = loc.Z();
121  out[3] = d1.X();
122  out[4] = d1.Y();
123  out[5] = d1.Z();
124  out[6] = d2.X();
125  out[7] = d2.Y();
126  out[8] = d2.Z();
127 
128  return out;
129 }
130 
131 Array<OneD, NekDouble> CADCurveOCE::NormalWRT(NekDouble t, int surf)
132 {
134  pair<CADSurfSharedPtr, CADOrientation::Orientation> surface;
135  ASSERTL0(m_adjSurfs.size() == 1, "This will only work in 2D for one surface at the moment");
136  surface = m_adjSurfs[0];
137 
138  Array<OneD, NekDouble> uv = surface.first->locuv(p);
139  Array<OneD, NekDouble> d1 = surface.first->D1(uv);
140 
141  NekDouble t1 = t - 1e-8;
142  NekDouble t2 = t + 1e-8;
143 
144  if(surface.second == CADOrientation::eBackwards)
145  {
146  swap(t1, t2);
147  }
148 
149  Array<OneD, NekDouble> uv1 = surface.first->locuv(P(t1));
150  Array<OneD, NekDouble> uv2 = surface.first->locuv(P(t2));
151 
152  NekDouble du = uv2[1] - uv1[1];
153  NekDouble dv = -1.0*(uv2[0] - uv1[0]);
154 
155  Array<OneD, NekDouble> N(3,0.0);
156  N[0] = (d1[3] * du + d1[6] * dv) / 2.0;
157  N[1] = (d1[4] * du + d1[7] * dv) / 2.0;
158  N[2] = (d1[5] * du + d1[8] * dv) / 2.0;
159 
160  NekDouble mag = sqrt(N[0]*N[0] + N[1]*N[1] + N[2]*N[2]);
161  N[0] /= mag;
162  N[1] /= mag;
163  N[2] /= mag;
164 
165  return N;
166 }
167 
169 {
170  GeomLProp_CLProps d(m_c,2,1e-8);
171  d.SetParameter(t+1e-8);
172 
173  gp_Vec d2 = d.D2();
174  if(d2.Magnitude() < 1e-8)
175  {
176  //no normal, stright line
177  return Array<OneD, NekDouble>(3,0.0);
178  }
179 
180  gp_Dir n;
181  d.Normal(n);
182 
184  N[0] = n.X();
185  N[1] = n.Y();
186  N[2] = n.Z();
187 
188  return N;
189 }
190 
191 NekDouble CADCurveOCE::Curvature(NekDouble t)
192 {
193  GeomLProp_CLProps d(m_c,2,1e-8);
194  d.SetParameter(t);
195 
196  return d.Curvature() * 1000.0;
197 }
198 
199 Array<OneD, NekDouble> CADCurveOCE::GetBounds()
200 {
202  t[0] = m_occCurve.FirstParameter();
203  t[1] = m_occCurve.LastParameter();
204 
205  return t;
206 }
207 
208 Array<OneD, NekDouble> CADCurveOCE::GetMinMax()
209 {
210  Array<OneD, NekDouble> locs(6);
211 
212  gp_Pnt start =
213  BRep_Tool::Pnt(TopExp::FirstVertex(m_occEdge, Standard_True));
214  gp_Pnt end = BRep_Tool::Pnt(TopExp::LastVertex(m_occEdge, Standard_True));
215 
216  locs[0] = start.X();
217  locs[1] = start.Y();
218  locs[2] = start.Z();
219  locs[3] = end.X();
220  locs[4] = end.Y();
221  locs[5] = end.Z();
222 
223  return locs;
224 }
225 
226 }
227 }
CADCurveFactory & GetCADCurveFactory()
Definition: CADSystem.cpp:64
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
STL namespace.
double NekDouble
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215