Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CurveMesh.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: CurveMesh.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: curvemesh object methods.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
38 
39 using namespace std;
40 namespace Nektar
41 {
42 namespace NekMeshUtils
43 {
44 
46 {
47  // this algorithm is mostly based on the work in chapter 19
48 
49  m_bounds = m_cadcurve->GetBounds();
50  m_curvelength = m_cadcurve->GetTotLength();
51  m_numSamplePoints =
52  int(m_curvelength / m_mesh->m_octree->GetMinDelta()) + 10;
53  ds = m_curvelength / (m_numSamplePoints - 1);
54 
55  GetSampleFunction();
56 
57  Ae = 0.0;
58 
59  for (int i = 0; i < m_numSamplePoints - 1; i++)
60  {
61  Ae += ds * (1.0 / m_dst[i][0] + 1.0 / m_dst[i + 1][0]) / 2.0;
62  }
63 
64  Ne = round(Ae);
65 
66  if (Ne + 1 < 2)
67  {
68  meshsvalue.resize(2);
69  meshsvalue[0] = 0.0;
70  meshsvalue[1] = m_curvelength;
71  Ne = 1;
72  }
73  else
74  {
75 
76  GetPhiFunction();
77 
78  meshsvalue.resize(Ne + 1);
79  meshsvalue[0] = 0.0;
80  meshsvalue[Ne] = m_curvelength;
81 
82  for (int i = 1; i < Ne; i++)
83  {
84  int iterationcounter = 0;
85  bool iterate = true;
86  int k = i;
87  NekDouble ski = meshsvalue[i - 1];
88  NekDouble lastSki;
89  while (iterate)
90  {
91  iterationcounter++;
92  NekDouble rhs = EvaluateDS(ski) / Ae * (EvaluatePS(ski) - k);
93  lastSki = ski;
94  ski = ski - rhs;
95  if (abs(lastSki - ski) < 1E-8)
96  {
97  iterate = false;
98  }
99 
100  ASSERTL0(iterationcounter < 1000000, "iteration failed");
101  }
102 
103  meshsvalue[i] = ski;
104  }
105  }
106 
107  NekDouble t;
109 
110  vector<CADVertSharedPtr> verts = m_cadcurve->GetVertex();
111  vector<pair<CADSurfSharedPtr, CADOrientation::Orientation> > s = m_cadcurve->GetAdjSurf();
112 
113  NodeSharedPtr n = verts[0]->GetNode();
114  t = m_bounds[0];
115  n->SetCADCurve(m_id, m_cadcurve, t);
116  loc = n->GetLoc();
117  for (int j = 0; j < s.size(); j++)
118  {
119  if (verts[0]->IsDegen() == s[j].first->GetId())
120  {
121  // if the degen has been set for this node the node
122  // already knows its corrected location
123  continue;
124  }
125 
126  Array<OneD, NekDouble> uv = s[j].first->locuv(loc);
127  n->SetCADSurf(s[j].first->GetId(), s[j].first, uv);
128  }
129  m_meshpoints.push_back(n);
130 
131  for (int i = 1; i < meshsvalue.size() - 1; i++)
132  {
133  t = m_cadcurve->tAtArcLength(meshsvalue[i]);
134  loc = m_cadcurve->P(t);
135  NodeSharedPtr n2 = boost::shared_ptr<Node>(
136  new Node(m_mesh->m_numNodes++, loc[0], loc[1], loc[2]));
137  n2->SetCADCurve(m_id, m_cadcurve, t);
138  for (int j = 0; j < s.size(); j++)
139  {
140  Array<OneD, NekDouble> uv = s[j].first->locuv(loc);
141  n2->SetCADSurf(s[j].first->GetId(), s[j].first, uv);
142  }
143  m_meshpoints.push_back(n2);
144  }
145 
146  n = verts[1]->GetNode();
147  t = m_bounds[1];
148  n->SetCADCurve(m_id, m_cadcurve, t);
149  loc = n->GetLoc();
150  for (int j = 0; j < s.size(); j++)
151  {
152  if (verts[1]->IsDegen() == s[j].first->GetId())
153  {
154  // if the degen has been set for this node the node
155  // already knows its corrected location
156  continue;
157  }
158 
159  Array<OneD, NekDouble> uv = s[j].first->locuv(loc);
160  n->SetCADSurf(s[j].first->GetId(), s[j].first, uv);
161  }
162  m_meshpoints.push_back(n);
163 
164  ASSERTL0(Ne + 1 == m_meshpoints.size(),
165  "incorrect number of points in curve mesh");
166 
167  // make edges and add them to the edgeset for the face mesher to use
168  for (int i = 0; i < m_meshpoints.size() - 1; i++)
169  {
170  EdgeSharedPtr e = boost::shared_ptr<Edge>(
171  new Edge(m_meshpoints[i], m_meshpoints[i + 1]));
172  e->m_parentCAD = m_cadcurve;
173  m_mesh->m_edgeSet.insert(e);
174  m_meshedges.push_back(e);
175  }
176 
177  if (m_mesh->m_verbose)
178  {
179  cout << "\r "
180  " "
181  " ";
182  cout << scientific << "\r\t\tCurve " << m_id << endl
183  << "\t\t\tLength: " << m_curvelength << endl
184  << "\t\t\tNodes: " << m_meshpoints.size() << endl
185  << "\t\t\tSample points: " << m_numSamplePoints << endl
186  << endl;
187  }
188 }
189 
190 void CurveMesh::GetPhiFunction()
191 {
192  m_ps.resize(m_numSamplePoints);
193  vector<NekDouble> newPhi;
194  newPhi.resize(2);
195 
196  newPhi[0] = 0.0;
197  newPhi[1] = 0.0;
198 
199  m_ps[0] = newPhi;
200 
201  NekDouble runningInt = 0.0;
202 
203  for (int i = 1; i < m_numSamplePoints; i++)
204  {
205  runningInt += (1.0 / m_dst[i - 1][0] + 1.0 / m_dst[i][0]) / 2.0 * ds;
206  newPhi[0] = Ne / Ae * runningInt;
207  newPhi[1] = m_dst[i][1];
208  m_ps[i] = newPhi;
209  }
210 }
211 
212 NekDouble CurveMesh::EvaluateDS(NekDouble s)
213 {
214  int a = 0;
215  int b = 0;
216 
217  ASSERTL1(!(s < 0) && !(s > m_curvelength),"s out of bounds");
218 
219  if (s == 0)
220  {
221  return m_dst[0][0];
222  }
223  else if (s == m_curvelength)
224  {
225  return m_dst[m_numSamplePoints - 1][0];
226  }
227 
228  for (int i = 0; i < m_numSamplePoints - 1; i++)
229  {
230  if (m_dst[i][1] < s && m_dst[i + 1][1] >= s)
231  {
232  a = i;
233  b = i + 1;
234  break;
235  }
236  }
237 
238  NekDouble s1 = m_dst[a][1];
239  NekDouble s2 = m_dst[b][1];
240  NekDouble d1 = m_dst[a][0];
241  NekDouble d2 = m_dst[b][0];
242 
243  NekDouble m = (d2 - d1) / (s2 - s1);
244  NekDouble c = d2 - m * s2;
245 
246  ASSERTL0(m * s + c == m * s + c, "DS"); // was getting nans here
247 
248  return m * s + c;
249 }
250 
251 NekDouble CurveMesh::EvaluatePS(NekDouble s)
252 {
253  int a = 0;
254  int b = 0;
255 
256  ASSERTL1(!(s < 0) && !(s > m_curvelength),"s out of bounds");
257 
258  if (s == 0)
259  {
260  return m_ps[0][0];
261  }
262  else if (s == m_curvelength)
263  {
264  return m_ps[m_numSamplePoints - 1][0];
265  }
266 
267  for (int i = 0; i < m_numSamplePoints - 1; i++)
268  {
269  if (m_ps[i][1] < s && m_ps[i + 1][1] >= s)
270  {
271  a = i;
272  b = i + 1;
273  break;
274  }
275  }
276 
277  if (a == b)
278  {
279  cout << endl;
280  cout << a << " " << b << endl;
281  cout << s << endl;
282  exit(-1);
283  }
284 
285  NekDouble s1 = m_ps[a][1];
286  NekDouble s2 = m_ps[b][1];
287  NekDouble d1 = m_ps[a][0];
288  NekDouble d2 = m_ps[b][0];
289 
290  NekDouble m = (d2 - d1) / (s2 - s1);
291  NekDouble c = d2 - m * s2;
292 
293  ASSERTL0(m * s + c == m * s + c, "PS");
294 
295  return m * s + c;
296 }
297 
298 void CurveMesh::GetSampleFunction()
299 {
300  m_dst.resize(m_numSamplePoints);
301 
302  vector<NekDouble> dsti;
303  dsti.resize(3);
304 
305  for (int i = 0; i < m_numSamplePoints; i++)
306  {
307  dsti[1] = i * ds;
308  NekDouble t = m_cadcurve->tAtArcLength(dsti[1]);
309 
310  Array<OneD, NekDouble> loc = m_cadcurve->P(t);
311 
312  /*NekDouble ts =
313  m_bl.Evaluate(m_blID, loc[0], loc[1], loc[2], 0.0);
314 
315  if (ts > 0.0)
316  {
317  Array<OneD, NekDouble> N = m_cadcurve->N(t);
318  Array<OneD, NekDouble> Nwrt = m_cadcurve->NormalWRT(t, 0);
319 
320  if(N[0]*N[0] + N[1]*N[1] + N[2]*N[2] < 1e-6)
321  {
322  dsti[0] = m_mesh->m_octree->Query(loc);
323  }
324  else if ( N[0]*Nwrt[0] + N[1]*Nwrt[1] + N[2]*Nwrt[2] > 0)
325  {
326  //concave
327  dsti[0] = m_mesh->m_octree->Query(loc);
328  }
329  else
330  {
331  NekDouble R = 1.0 / m_cadcurve->Curvature(t);
332  if(R > 2.0*t)
333  {
334  R = 2.0*t;
335  }
336  Array<OneD, NekDouble> tloc(3);
337  tloc[0] = loc[0] + ts * Nwrt[0];
338  tloc[1] = loc[1] + ts * Nwrt[1];
339  tloc[2] = loc[2] + ts * Nwrt[2];
340 
341  NekDouble d = m_mesh->m_octree->Query(tloc);
342 
343  dsti[0] = d * R / (R + ts);
344  }
345  }
346  else
347  {*/
348  dsti[0] = m_mesh->m_octree->Query(loc);
349  //}
350 
351  dsti[2] = t;
352 
353  m_dst[i] = dsti;
354  }
355 }
356 }
357 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Represents an edge which joins two points.
Definition: Edge.h:58
STL namespace.
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
double NekDouble
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:136
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228