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