Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
FaceMesh.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: SurfaceMesh.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: surfacemesh object methods.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
39 #include <limits>
40 
41 using namespace std;
42 namespace Nektar
43 {
44 namespace NekMeshUtils
45 {
46 
47 bool FaceMesh::ValidateCurves()
48 {
49  vector<int> curvesInSurface;
50  for (int i = 0; i < m_edgeloops.size(); i++)
51  {
52  for (int j = 0; j < m_edgeloops[i]->edges.size(); j++)
53  {
54  curvesInSurface.push_back(m_edgeloops[i]->edges[j]->GetId());
55  }
56  }
57 
58  bool error = false;
59 
60  for (int i = 0; i < curvesInSurface.size(); i++)
61  {
62  vector<EdgeSharedPtr> es =
63  m_curvemeshes[curvesInSurface[i]]->GetMeshEdges();
64 
65  for (int j = i; j < curvesInSurface.size(); j++)
66  {
67  if (i == j)
68  {
69  continue;
70  }
71 
72  vector<EdgeSharedPtr> es2 =
73  m_curvemeshes[curvesInSurface[j]]->GetMeshEdges();
74 
75  for (int l = 0; l < es.size(); l++)
76  {
77  Array<OneD, NekDouble> P1 = es[l]->m_n1->GetCADSurfInfo(m_id);
78  Array<OneD, NekDouble> P2 = es[l]->m_n2->GetCADSurfInfo(m_id);
79  for (int k = 0; k < es2.size(); k++)
80  {
81  if (es[l]->m_n1 == es2[k]->m_n1 ||
82  es[l]->m_n1 == es2[k]->m_n2 ||
83  es[l]->m_n2 == es2[k]->m_n1 ||
84  es[l]->m_n2 == es2[k]->m_n2)
85  {
86  continue;
87  }
88 
90  es2[k]->m_n1->GetCADSurfInfo(m_id);
92  es2[k]->m_n2->GetCADSurfInfo(m_id);
93 
94  NekDouble den = (P4[0] - P3[0]) * (P2[1] - P1[1]) -
95  (P2[0] - P1[0]) * (P4[1] - P3[1]);
96 
97  if (fabs(den) < 1e-8)
98  {
99  continue;
100  }
101 
102  NekDouble t = ((P1[0] - P3[0]) * (P4[1] - P3[1]) -
103  (P4[0] - P3[0]) * (P1[1] - P3[1])) /
104  den;
105  NekDouble u =
106  (P1[0] - P3[0] + t * (P2[0] - P1[0])) / (P4[0] - P3[0]);
107 
108  if (t < 1.0 && t > 0.0 && u < 1.0 && u > 0.0)
109  {
111  uv[0] = P1[0] + t * (P2[0] - P1[0]);
112  uv[1] = P1[1] + t * (P2[1] - P1[1]);
113  Array<OneD, NekDouble> loc = m_cadsurf->P(uv);
114  cout << endl
115  << "Curve mesh error at " << loc[0] << " "
116  << loc[1] << " " << loc[2] << " on face " << m_id
117  << endl;
118  error = true;
119  }
120  }
121  }
122  }
123  }
124  return error;
125 }
126 
128 {
129  Stretching();
130  OrientateCurves();
131 
132  int numPoints = 0;
133  for (int i = 0; i < orderedLoops.size(); i++)
134  {
135  numPoints += orderedLoops[i].size();
136  for (int j = 0; j < orderedLoops[i].size(); j++)
137  {
138  m_inBoundary.insert(orderedLoops[i][j]);
139  }
140  }
141 
142  stringstream ss;
143  ss << "3 points required for triangulation, " << numPoints << " in loop"
144  << endl;
145  ss << "curves: ";
146  for (int i = 0; i < m_edgeloops.size(); i++)
147  {
148  for (int j = 0; j < m_edgeloops[i]->edges.size(); j++)
149  {
150  ss << m_edgeloops[i]->edges[j]->GetId() << " ";
151  }
152  }
153 
154  ASSERTL0(numPoints > 2, ss.str());
155 
156  // create interface to triangle thirdparty library
157  TriangleInterfaceSharedPtr pplanemesh =
159 
160  vector<Array<OneD, NekDouble> > centers;
161  for (int i = 0; i < m_edgeloops.size(); i++)
162  {
163  centers.push_back(m_edgeloops[i]->center);
164  }
165 
166  pplanemesh->Assign(orderedLoops, centers, m_id, m_str);
167 
168  pplanemesh->Mesh();
169 
170  pplanemesh->Extract(m_connec);
171 
172  bool repeat = true;
173  int meshcounter = 1;
174 
175  // continuously remesh until all triangles conform to the spacing in the
176  // octree
177  while (repeat)
178  {
179  repeat = Validate();
180  if (!repeat)
181  {
182  break;
183  }
184  m_connec.clear();
185  pplanemesh->AssignStiener(m_stienerpoints);
186  pplanemesh->Mesh();
187  pplanemesh->Extract(m_connec);
188  meshcounter++;
189  // break;
190  }
191 
192  // build a local version of the mesh (one set of triangles). this is done
193  // so edge connectivity infomration can be used for optimisation
194  BuildLocalMesh();
195 
196  OptimiseLocalMesh();
197 
198  // make new elements and add to list from list of nodes and connectivity
199  // from triangle removing unnesercary infomration from the elements
200  for (int i = 0; i < m_localElements.size(); i++)
201  {
202  vector<EdgeSharedPtr> e = m_localElements[i]->GetEdgeList();
203  for (int j = 0; j < e.size(); j++)
204  {
205  e[j]->m_elLink.clear();
206  }
207  m_mesh->m_element[2].push_back(m_localElements[i]);
208  }
209 
210  if (m_mesh->m_verbose)
211  {
212  cout << "\r "
213  " "
214  " ";
215  cout << scientific << "\r\t\tFace " << m_id << endl
216  << "\t\t\tNodes: " << m_localNodes.size() << endl
217  << "\t\t\tEdges: " << m_localEdges.size() << endl
218  << "\t\t\tTriangles: " << m_localElements.size() << endl
219  << "\t\t\tLoops: " << m_edgeloops.size() << endl
220  << "\t\t\tSTR: " << m_str << endl
221  << endl;
222  }
223 }
224 
225 void FaceMesh::OptimiseLocalMesh()
226 {
227  // each optimisation algorithm is based on the work in chapter 19
228  DiagonalSwap();
229 
230  Smoothing();
231 
232  DiagonalSwap();
233 
234  Smoothing();
235 }
236 
237 void FaceMesh::Smoothing()
238 {
239  EdgeSet::iterator eit;
240  NodeSet::iterator nit;
241 
242  Array<OneD, NekDouble> bounds = m_cadsurf->GetBounds();
243 
244  map<int, vector<EdgeSharedPtr> > connectingedges;
245 
246  map<int, vector<ElementSharedPtr> > connectingelements;
247 
248  for (eit = m_localEdges.begin(); eit != m_localEdges.end(); eit++)
249  {
250  connectingedges[(*eit)->m_n1->m_id].push_back(*eit);
251  connectingedges[(*eit)->m_n2->m_id].push_back(*eit);
252  }
253 
254  for (int i = 0; i < m_localElements.size(); i++)
255  {
256  vector<NodeSharedPtr> v = m_localElements[i]->GetVertexList();
257  for (int j = 0; j < 3; j++)
258  {
259  connectingelements[v[j]->m_id].push_back(m_localElements[i]);
260  }
261  }
262 
263  // perform 4 runs of elastic relaxation based on the octree
264  for (int q = 0; q < 4; q++)
265  {
266  for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
267  {
268  NodeSet::iterator f = m_inBoundary.find((*nit));
269 
270  // node is on curve so skip
271  if (f != m_inBoundary.end())
272  {
273  continue;
274  }
275 
276  // this can be real nodes or dummy nodes depending on the system
277  vector<NodeSharedPtr> connodes;
278 
279  vector<EdgeSharedPtr> edges = connectingedges[(*nit)->m_id];
280  vector<ElementSharedPtr> els = connectingelements[(*nit)->m_id];
281 
282  vector<NodeSharedPtr> nodesystem;
283  vector<NekDouble> lamp;
284 
285  for (int i = 0; i < edges.size(); i++)
286  {
287  vector<NekDouble> lambda;
288 
289  NodeSharedPtr J;
290  if (*nit == edges[i]->m_n1)
291  {
292  J = edges[i]->m_n2;
293  }
294  else if (*nit == edges[i]->m_n2)
295  {
296  J = edges[i]->m_n1;
297  }
298  else
299  {
300  ASSERTL0(false, "could not find node");
301  }
302 
303  Array<OneD, NekDouble> ui = (*nit)->GetCADSurfInfo(m_id);
304  Array<OneD, NekDouble> uj = J->GetCADSurfInfo(m_id);
305 
306  for (int j = 0; j < els.size(); j++)
307  {
308  vector<NodeSharedPtr> v = els[j]->GetVertexList();
309 
310  // elememt is adjacent to J therefore no intersection on IJ
311  if (v[0] == J || v[1] == J || v[2] == J)
312  {
313  continue;
314  }
315 
316  // need to find other edge
317  EdgeSharedPtr AtoB;
318  bool found = false;
319  vector<EdgeSharedPtr> es = els[j]->GetEdgeList();
320  for (int k = 0; k < es.size(); k++)
321  {
322  if (!(es[k]->m_n1 == *nit || es[k]->m_n2 == *nit))
323  {
324  found = true;
325  AtoB = es[k];
326  break;
327  }
328  }
329  ASSERTL0(found, "failed to find edge to test");
330 
331  Array<OneD, NekDouble> A = AtoB->m_n1->GetCADSurfInfo(m_id);
332  Array<OneD, NekDouble> B = AtoB->m_n2->GetCADSurfInfo(m_id);
333 
334  NekDouble lam = ((A[0] - uj[0]) * (B[1] - A[1]) -
335  (A[1] - uj[1]) * (B[0] - A[0])) /
336  ((ui[0] - uj[0]) * (B[1] - A[1]) -
337  (ui[1] - uj[1]) * (B[0] - A[0]));
338 
339  if (!(lam < 0) && !(lam > 1))
340  lambda.push_back(lam);
341  }
342 
343  if (lambda.size() > 0)
344  {
345  sort(lambda.begin(), lambda.end());
346  // make a new dummy node based on the system
348  ud[0] = uj[0] + lambda[0] * (ui[0] - uj[0]);
349  ud[1] = uj[1] + lambda[0] * (ui[1] - uj[1]);
350  Array<OneD, NekDouble> locd = m_cadsurf->P(ud);
351  NodeSharedPtr dn = boost::shared_ptr<Node>(
352  new Node(0, locd[0], locd[1], locd[2]));
353  dn->SetCADSurf(m_id, m_cadsurf, ud);
354 
355  nodesystem.push_back(dn);
356  lamp.push_back(lambda[0]);
357  }
358  else
359  {
360  nodesystem.push_back(J);
361  lamp.push_back(1.0);
362  }
363  }
364 
366  u0[0] = 0.0;
367  u0[1] = 0.0;
368 
369  for (int i = 0; i < nodesystem.size(); i++)
370  {
371  Array<OneD, NekDouble> uj = nodesystem[i]->GetCADSurfInfo(m_id);
372  u0[0] += uj[0] / nodesystem.size();
373  u0[1] += uj[1] / nodesystem.size();
374  }
375 
376  /*Array<OneD, NekDouble> pu0 = m_cadsurf->P(u0);
377  NekDouble di = m_mesh->m_octree->Query(pu0);
378  Array<OneD, NekDouble> F(2, 0.0), dF(4, 0.0);
379  for (int i = 0; i < nodesystem.size(); i++)
380  {
381  Array<OneD, NekDouble> rj = nodesystem[i]->GetLoc();
382  Array<OneD, NekDouble> uj = nodesystem[i]->GetCADSurfInfo(m_id);
383  NekDouble dj = m_mesh->m_octree->Query(rj);
384  NekDouble d = (di + dj) / 2.0;
385  NekDouble wij = sqrt((rj[0] - pu0[0]) * (rj[0] - pu0[0]) +
386  (rj[1] - pu0[1]) * (rj[1] - pu0[1]) +
387  (rj[2] - pu0[2]) * (rj[2] - pu0[2])) -
388  d;
389 
390  NekDouble umag = sqrt((uj[0] - u0[0]) * (uj[0] - u0[0]) +
391  (uj[1] - u0[1]) * (uj[1] - u0[1]));
392 
393  F[0] += wij * (uj[0] - u0[0]) / umag;
394  F[1] += wij * (uj[1] - u0[1]) / umag;
395 
396  Array<OneD, NekDouble> d1 = m_cadsurf->D1(u0);
397 
398  Array<OneD, NekDouble> dw(2, 0.0);
399  dw[0] = -2.0 *
400  ((rj[0] - pu0[0]) * d1[3] + (rj[1] - pu0[1]) * d1[4] +
401  (rj[2] - pu0[2]) * d1[5]);
402  dw[1] = -2.0 *
403  ((rj[0] - pu0[0]) * d1[6] + (rj[1] - pu0[1]) * d1[7] +
404  (rj[2] - pu0[2]) * d1[8]);
405 
406  Array<OneD, NekDouble> drhs(4, 0.0);
407  drhs[0] = 2.0 * ((uj[0] - u0[0]) * (uj[0] - u0[0]) - umag) /
408  umag / umag;
409  drhs[1] =
410  2.0 * ((uj[0] - u0[0]) * (uj[1] - u0[1])) / umag / umag;
411  drhs[2] = drhs[1];
412  drhs[3] = 2.0 * ((uj[1] - u0[1]) * (uj[1] - u0[1]) - umag) /
413  umag / umag;
414 
415  dF[0] += dw[0] * (uj[0] - u0[0]) / umag + wij * drhs[0];
416 
417  dF[1] += dw[0] * (uj[1] - u0[1]) / umag + wij * drhs[1];
418 
419  dF[2] += dw[1] * (uj[0] - u0[0]) / umag + wij * drhs[2];
420 
421  dF[3] += dw[1] * (uj[1] - u0[1]) / umag + wij * drhs[3];
422  }
423 
424  NekDouble det = dF[0] * dF[3] - dF[1] * dF[2];
425  NekDouble tmp = dF[0];
426  dF[0] = dF[3] / det;
427  dF[3] = tmp / det;
428  dF[1] *= -1.0 / det;
429  dF[2] *= -1.0 / det;
430 
431  u0[0] -= (dF[0] * F[0] + dF[2] * F[1]);
432  u0[1] -= (dF[1] * F[0] + dF[3] * F[1]);*/
433 
434  bool inbounds = true;
435  if (u0[0] < bounds[0])
436  {
437  inbounds = false;
438  }
439  else if (u0[0] > bounds[1])
440  {
441  inbounds = false;
442  }
443  else if (u0[1] < bounds[2])
444  {
445  inbounds = false;
446  }
447  else if (u0[1] > bounds[3])
448  {
449  inbounds = false;
450  }
451 
452  if (!inbounds)
453  {
454  continue;
455  }
456 
457  /*Array<OneD, NekDouble> FN(2, 0.0);
458  pu0 = m_cadsurf->P(u0);
459  di = m_mesh->m_octree->Query(pu0);
460  for (int i = 0; i < nodesystem.size(); i++)
461  {
462  Array<OneD, NekDouble> rj = nodesystem[i]->GetLoc();
463  Array<OneD, NekDouble> uj = nodesystem[i]->GetCADSurfInfo(m_id);
464  NekDouble dj = m_mesh->m_octree->Query(rj);
465  NekDouble d = (di + dj) / 2.0;
466  NekDouble wij = sqrt((rj[0] - pu0[0]) * (rj[0] - pu0[0]) +
467  (rj[1] - pu0[1]) * (rj[1] - pu0[1]) +
468  (rj[2] - pu0[2]) * (rj[2] - pu0[2])) - d;
469 
470  NekDouble umag = sqrt((uj[0] - u0[0]) * (uj[0] - u0[0]) +
471  (uj[1] - u0[1]) * (uj[1] - u0[1]));
472 
473  FN[0] += wij * (uj[0] - u0[0]) / umag;
474  FN[1] += wij * (uj[1] - u0[1]) / umag;
475  }
476 
477  if (F[0] * F[0] + F[1] * F[1] < FN[0] * FN[0] + FN[1] * FN[1])
478  {
479  continue;
480  }*/
481 
482  Array<OneD, NekDouble> l2 = m_cadsurf->P(u0);
483  (*nit)->Move(l2, m_id, u0);
484  }
485  }
486 }
487 
488 void FaceMesh::DiagonalSwap()
489 {
490  map<int, int> idealConnec;
491  map<int, int> actualConnec;
492  map<int, vector<EdgeSharedPtr> > nodetoedge;
493  // figure out ideal node count and actual node count
494  EdgeSet::iterator eit;
495  for (eit = m_localEdges.begin(); eit != m_localEdges.end(); eit++)
496  {
497  nodetoedge[(*eit)->m_n1->m_id].push_back(*eit);
498  nodetoedge[(*eit)->m_n2->m_id].push_back(*eit);
499  }
500  NodeSet::iterator nit;
501  for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
502  {
503  // this routine is broken and needs looking at
504  NodeSet::iterator f = m_inBoundary.find((*nit));
505  if (f == m_inBoundary.end()) // node isnt on curve so skip
506  {
507  // node is interior
508  idealConnec[(*nit)->m_id] = 6;
509  }
510  else
511  {
512  // need to identify the two other nodes on the boundary to find
513  // interior angle
514  /*vector<NodeSharedPtr> ns;
515  vector<EdgeSharedPtr> e = nodetoedge[(*nit)->m_id];
516  for (int i = 0; i < e.size(); i++)
517  {
518  if (!e[i]->onCurve)
519  continue; // the linking nodes are not going to exist on
520  // interior edges
521 
522  if (e[i]->m_n1 == (*nit))
523  ns.push_back(e[i]->m_n2);
524  else
525  ns.push_back(e[i]->m_n1);
526  }
527  ASSERTL0(ns.size() == 2,
528  "failed to find 2 nodes in the angle system");
529 
530  idealConnec[(*nit)->m_id] =
531  ceil((*nit)->Angle(ns[0], ns[1]) / 3.142 * 3) + 1;*/
532  idealConnec[(*nit)->m_id] = 4;
533  }
534  }
535  for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
536  {
537  actualConnec[(*nit)->m_id] = nodetoedge[(*nit)->m_id].size();
538  }
539 
540  // edgeswapping fun times
541  // perfrom edge swap based on node defect and then angle
542  for (int q = 0; q < 4; q++)
543  {
544  int edgesStart = m_localEdges.size();
545  EdgeSet edges = m_localEdges;
546  m_localEdges.clear();
547 
548  int swappedEdges = 0;
549 
551 
552  for (it = edges.begin(); it != edges.end(); it++)
553  {
554  EdgeSharedPtr e = *it;
555 
556  NodeSet::iterator f1 = m_inBoundary.find((*it)->m_n1);
557  NodeSet::iterator f2 = m_inBoundary.find((*it)->m_n2);
558  if (f1 != m_inBoundary.end() && f2 != m_inBoundary.end())
559  {
560  m_localEdges.insert(e);
561  continue;
562  }
563 
564  ElementSharedPtr tri1 = e->m_elLink[0].first;
565  ElementSharedPtr tri2 = e->m_elLink[1].first;
566 
567  NodeSharedPtr n1 = e->m_n1;
568  NodeSharedPtr n2 = e->m_n2;
569 
570  vector<NodeSharedPtr> nt = tri1->GetVertexList();
571 
572  // identify node a,b,c,d of the swapping
573  NodeSharedPtr A, B, C, D;
574  if (nt[0] != n1 && nt[0] != n2)
575  {
576  C = nt[0];
577  B = nt[1];
578  A = nt[2];
579  }
580  else if (nt[1] != n1 && nt[1] != n2)
581  {
582  C = nt[1];
583  B = nt[2];
584  A = nt[0];
585  }
586  else if (nt[2] != n1 && nt[2] != n2)
587  {
588  C = nt[2];
589  B = nt[0];
590  A = nt[1];
591  }
592  else
593  {
594  ASSERTL0(false, "failed to identify verticies in tri1");
595  }
596 
597  nt = tri2->GetVertexList();
598 
599  if (nt[0] != n1 && nt[0] != n2)
600  {
601  D = nt[0];
602  }
603  else if (nt[1] != n1 && nt[1] != n2)
604  {
605  D = nt[1];
606  }
607  else if (nt[2] != n1 && nt[2] != n2)
608  {
609  D = nt[2];
610  }
611  else
612  {
613  ASSERTL0(false, "failed to identify verticies in tri2");
614  }
615 
616  // determine signed area of alternate config
617  // //cout<< A->GetNumCADSurf() << " " << B->GetNumCADSurf() << " "
618  // // << C->GetNumCADSurf() << " " << D->GetNumCADSurf() << endl;
619  // ofstream file;
620  // file.open("pts.3D");
621  // file << "x y z value" << endl;
622  // file << A->m_x << " " << A->m_y << " " << A->m_z << endl;
623  // file << B->m_x << " " << B->m_y << " " << B->m_z << endl;
624  // file << C->m_x << " " << C->m_y << " " << C->m_z << endl;
625  // file << D->m_x << " " << D->m_y << " " << D->m_z << endl;
626  // file.close();
627  Array<OneD, NekDouble> ai, bi, ci, di;
628  ai = A->GetCADSurfInfo(m_id);
629  bi = B->GetCADSurfInfo(m_id);
630  ci = C->GetCADSurfInfo(m_id);
631  di = D->GetCADSurfInfo(m_id);
632 
633  NekDouble CDA, CBD;
634 
635  CDA = 0.5 * (-di[0] * ci[1] + ai[0] * ci[1] + ci[0] * di[1] -
636  ai[0] * di[1] - ci[0] * ai[1] + di[0] * ai[1]);
637 
638  CBD = 0.5 * (-bi[0] * ci[1] + di[0] * ci[1] + ci[0] * bi[1] -
639  di[0] * bi[1] - ci[0] * di[1] + bi[0] * di[1]);
640 
641  // if signed area of the swapping triangles is less than zero
642  // that configuration is invalid and swap cannot be performed
643  if (!(CDA > 0.001 && CBD > 0.001))
644  {
645  m_localEdges.insert(e);
646  continue;
647  }
648 
649  bool swap = false; // assume do not swap
650 
651  if (q < 2)
652  {
653  int nodedefectbefore = 0;
654  nodedefectbefore +=
655  abs(actualConnec[A->m_id] - idealConnec[A->m_id]);
656  nodedefectbefore +=
657  abs(actualConnec[B->m_id] - idealConnec[B->m_id]);
658  nodedefectbefore +=
659  abs(actualConnec[C->m_id] - idealConnec[C->m_id]);
660  nodedefectbefore +=
661  abs(actualConnec[D->m_id] - idealConnec[D->m_id]);
662 
663  int nodedefectafter = 0;
664  nodedefectafter +=
665  abs(actualConnec[A->m_id] - 1 - idealConnec[A->m_id]);
666  nodedefectafter +=
667  abs(actualConnec[B->m_id] - 1 - idealConnec[B->m_id]);
668  nodedefectafter +=
669  abs(actualConnec[C->m_id] + 1 - idealConnec[C->m_id]);
670  nodedefectafter +=
671  abs(actualConnec[D->m_id] + 1 - idealConnec[D->m_id]);
672 
673  if (nodedefectafter < nodedefectbefore)
674  {
675  swap = true;
676  }
677  }
678  else
679  {
680  NekDouble minanglebefore = C->Angle(A, B);
681  minanglebefore = min(minanglebefore, A->Angle(B, C));
682  minanglebefore = min(minanglebefore, B->Angle(A, C));
683  minanglebefore = min(minanglebefore, B->Angle(A, D));
684  minanglebefore = min(minanglebefore, A->Angle(B, D));
685  minanglebefore = min(minanglebefore, D->Angle(A, B));
686 
687  NekDouble minangleafter = C->Angle(B, D);
688  minangleafter = min(minangleafter, D->Angle(B, C));
689  minangleafter = min(minangleafter, B->Angle(C, D));
690  minangleafter = min(minangleafter, C->Angle(A, D));
691  minangleafter = min(minangleafter, A->Angle(C, D));
692  minangleafter = min(minangleafter, D->Angle(A, C));
693 
694  if (minangleafter > minanglebefore)
695  {
696  swap = true;
697  }
698  }
699 
700  if (swap)
701  {
702  actualConnec[A->m_id]--;
703  actualConnec[B->m_id]--;
704  actualConnec[C->m_id]++;
705  actualConnec[D->m_id]++;
706 
707  // make the 4 other edges
708  EdgeSharedPtr CA, AD, DB, BC, CAt, ADt, DBt, BCt;
709  CAt = boost::shared_ptr<Edge>(new Edge(C, A));
710  ADt = boost::shared_ptr<Edge>(new Edge(A, D));
711  DBt = boost::shared_ptr<Edge>(new Edge(D, B));
712  BCt = boost::shared_ptr<Edge>(new Edge(B, C));
713 
714  vector<EdgeSharedPtr> es = tri1->GetEdgeList();
715  for (int i = 0; i < 3; i++)
716  {
717  if (es[i] == CAt)
718  {
719  CA = es[i];
720  }
721  if (es[i] == BCt)
722  {
723  BC = es[i];
724  }
725  }
726  es = tri2->GetEdgeList();
727  for (int i = 0; i < 3; i++)
728  {
729  if (es[i] == DBt)
730  {
731  DB = es[i];
732  }
733  if (es[i] == ADt)
734  {
735  AD = es[i];
736  }
737  }
738 
739  // now sort out links for the 4 edges surrounding the patch
740  vector<pair<ElementSharedPtr, int> > links;
741 
742  links = CA->m_elLink;
743  CA->m_elLink.clear();
744  for (int i = 0; i < links.size(); i++)
745  {
746  if (links[i].first->GetId() == tri1->GetId())
747  {
748  continue;
749  }
750  CA->m_elLink.push_back(links[i]);
751  }
752 
753  links = BC->m_elLink;
754  BC->m_elLink.clear();
755  for (int i = 0; i < links.size(); i++)
756  {
757  if (links[i].first->GetId() == tri1->GetId())
758  {
759  continue;
760  }
761  BC->m_elLink.push_back(links[i]);
762  }
763 
764  links = AD->m_elLink;
765  AD->m_elLink.clear();
766  for (int i = 0; i < links.size(); i++)
767  {
768  if (links[i].first->GetId() == tri2->GetId())
769  {
770  continue;
771  }
772  AD->m_elLink.push_back(links[i]);
773  }
774 
775  links = DB->m_elLink;
776  DB->m_elLink.clear();
777  for (int i = 0; i < links.size(); i++)
778  {
779  if (links[i].first->GetId() == tri2->GetId())
780  {
781  continue;
782  }
783  DB->m_elLink.push_back(links[i]);
784  }
785 
786  EdgeSharedPtr newe = boost::shared_ptr<Edge>(new Edge(C, D));
787 
788  vector<NodeSharedPtr> t1, t2;
789  t1.push_back(B);
790  t1.push_back(D);
791  t1.push_back(C);
792  t2.push_back(A);
793  t2.push_back(C);
794  t2.push_back(D);
795 
796  ElmtConfig conf(LibUtilities::eTriangle, 1, false, false);
797  vector<int> tags = tri1->GetTagList();
798 
799  int id1 = tri1->GetId();
800  int id2 = tri2->GetId();
801 
803  LibUtilities::eTriangle, conf, t1, tags);
804  tags = tri2->GetTagList();
806  LibUtilities::eTriangle, conf, t2, tags);
807 
808  ntri1->SetId(id1);
809  ntri2->SetId(id2);
810  ntri1->m_parentCAD = m_cadsurf;
811  ntri2->m_parentCAD = m_cadsurf;
812 
813  vector<EdgeSharedPtr> t1es = ntri1->GetEdgeList();
814  for (int i = 0; i < 3; i++)
815  {
816  if (t1es[i] == DB)
817  {
818  ntri1->SetEdge(i, DB);
819  DB->m_elLink.push_back(
820  pair<ElementSharedPtr, int>(ntri1, i));
821  }
822  else if (t1es[i] == BC)
823  {
824  ntri1->SetEdge(i, BC);
825  BC->m_elLink.push_back(
826  pair<ElementSharedPtr, int>(ntri1, i));
827  }
828  else if (t1es[i] == newe)
829  {
830  ntri1->SetEdge(i, newe);
831  newe->m_elLink.push_back(
832  pair<ElementSharedPtr, int>(ntri1, i));
833  }
834  else
835  {
836  ASSERTL0(false, "weird edge in new tri 1");
837  }
838  }
839  vector<EdgeSharedPtr> t2es = ntri2->GetEdgeList();
840  for (int i = 0; i < 3; i++)
841  {
842  if (t2es[i] == CA)
843  {
844  ntri2->SetEdge(i, CA);
845  CA->m_elLink.push_back(
846  pair<ElementSharedPtr, int>(ntri2, i));
847  }
848  else if (t2es[i] == AD)
849  {
850  ntri2->SetEdge(i, AD);
851  AD->m_elLink.push_back(
852  pair<ElementSharedPtr, int>(ntri2, i));
853  }
854  else if (t2es[i] == newe)
855  {
856  ntri2->SetEdge(i, newe);
857  newe->m_elLink.push_back(
858  pair<ElementSharedPtr, int>(ntri2, i));
859  }
860  else
861  {
862  ASSERTL0(false, "weird edge in new tri 2");
863  }
864  }
865 
866  m_localEdges.insert(newe);
867 
868  m_localElements[id1] = ntri1;
869  m_localElements[id2] = ntri2;
870 
871  swappedEdges++;
872  }
873  else
874  {
875  m_localEdges.insert(e);
876  }
877  }
878 
879  ASSERTL0(m_localEdges.size() == edgesStart, "mismatch edge count");
880  }
881 }
882 
883 void FaceMesh::BuildLocalMesh()
884 {
885  /*************************
886  // build a local set of nodes edges and elemenets for optimstaion prior to
887  putting them into m_mesh
888  */
889 
890  for (int i = 0; i < m_connec.size(); i++)
891  {
892  ElmtConfig conf(LibUtilities::eTriangle, 1, false, false);
893 
894  vector<int> tags;
895  tags.push_back(m_compId);
897  LibUtilities::eTriangle, conf, m_connec[i], tags);
898  E->m_parentCAD = m_cadsurf;
899 
900  vector<NodeSharedPtr> nods = E->GetVertexList();
901  for (int j = 0; j < nods.size(); j++)
902  {
903  // nodes are already unique some will insert some wont
904  m_localNodes.insert(nods[j]);
905  }
906 
907  E->SetId(m_localElements.size());
908  m_localElements.push_back(E);
909  }
910 
911  for (int i = 0; i < m_localElements.size(); ++i)
912  {
913  for (int j = 0; j < m_localElements[i]->GetEdgeCount(); ++j)
914  {
915  pair<EdgeSet::iterator, bool> testIns;
916  EdgeSharedPtr ed = m_localElements[i]->GetEdge(j);
917  // look for edge in m_mesh edgeset from curves
918  EdgeSet::iterator s = m_mesh->m_edgeSet.find(ed);
919  if (!(s == m_mesh->m_edgeSet.end()))
920  {
921  ed = *s;
922  m_localElements[i]->SetEdge(j, *s);
923  }
924 
925  testIns = m_localEdges.insert(ed);
926 
927  if (testIns.second)
928  {
929  EdgeSharedPtr ed2 = *testIns.first;
930  ed2->m_elLink.push_back(
931  pair<ElementSharedPtr, int>(m_localElements[i], j));
932  }
933  else
934  {
935  EdgeSharedPtr e2 = *(testIns.first);
936  m_localElements[i]->SetEdge(j, e2);
937  e2->m_elLink.push_back(
938  pair<ElementSharedPtr, int>(m_localElements[i], j));
939  }
940  }
941  }
942 }
943 
944 void FaceMesh::Stretching()
945 {
946  // define a sampling and calculate the aspect ratio of the paramter plane
947  m_str = 0.0;
948  Array<OneD, NekDouble> bnds = m_cadsurf->GetBounds();
949 
950  NekDouble dxu = int(bnds[1] - bnds[0] < bnds[3] - bnds[2]
951  ? 40
952  : (bnds[1] - bnds[0]) / (bnds[3] - bnds[2]) * 40);
953  NekDouble dxv = int(bnds[3] - bnds[2] < bnds[1] - bnds[0]
954  ? 40
955  : (bnds[3] - bnds[2]) / (bnds[1] - bnds[0]) * 40);
956 
957  NekDouble du = (bnds[1] - bnds[0]) / dxu;
958  NekDouble dv = (bnds[3] - bnds[2]) / dxv;
959 
960  int ct = 0;
961 
962  for (int i = 0; i < dxu; i++)
963  {
964  for (int j = 0; j < dxv; j++)
965  {
967  uv[0] = bnds[0] + i * du;
968  uv[1] = bnds[2] + j * dv;
969  if (i == dxu - 1)
970  {
971  uv[0] = bnds[1];
972  }
973  if (j == dxv - 1)
974  {
975  uv[1] = bnds[3];
976  }
977  Array<OneD, NekDouble> r = m_cadsurf->D1(uv);
978 
979  NekDouble ru = sqrt(r[3] * r[3] + r[4] * r[4] + r[5] * r[5]);
980  NekDouble rv = sqrt(r[6] * r[6] + r[7] * r[7] + r[8] * r[8]);
981 
982  ru *= du;
983  rv *= dv;
984 
985  if (rv < 1E-8)
986  {
987  continue;
988  }
989 
990  m_str += ru / rv;
991  ct++;
992  }
993  }
994 
995  m_str /= ct;
996 }
997 
998 bool FaceMesh::Validate()
999 {
1000  // check all edges in the current mesh for length against the octree
1001  // if the octree is not conformed to add a new point inside the triangle
1002  // if no new points are added meshing can stop
1003  int pointBefore = m_stienerpoints.size();
1004  for (int i = 0; i < m_connec.size(); i++)
1005  {
1006  Array<OneD, NekDouble> r(3), a(3);
1007 
1008  vector<Array<OneD, NekDouble> > info;
1009 
1010  for (int j = 0; j < 3; j++)
1011  {
1012  info.push_back(m_connec[i][j]->GetCADSurfInfo(m_id));
1013  }
1014 
1015  r[0] = m_connec[i][0]->Distance(m_connec[i][1]);
1016  r[1] = m_connec[i][1]->Distance(m_connec[i][2]);
1017  r[2] = m_connec[i][2]->Distance(m_connec[i][0]);
1018 
1019  a[0] = m_connec[i][0]->Angle(m_connec[i][1]->GetLoc(),
1020  m_connec[i][2]->GetLoc(),
1021  m_cadsurf->N(info[0]));
1022  a[1] = m_connec[i][1]->Angle(m_connec[i][2]->GetLoc(),
1023  m_connec[i][0]->GetLoc(),
1024  m_cadsurf->N(info[1]));
1025  a[2] = m_connec[i][2]->Angle(m_connec[i][0]->GetLoc(),
1026  m_connec[i][1]->GetLoc(),
1027  m_cadsurf->N(info[2]));
1028 
1029  NekDouble d1 = m_mesh->m_octree->Query(m_connec[i][0]->GetLoc());
1030  NekDouble d2 = m_mesh->m_octree->Query(m_connec[i][1]->GetLoc());
1031  NekDouble d3 = m_mesh->m_octree->Query(m_connec[i][2]->GetLoc());
1032 
1033  Array<OneD, NekDouble> uvc(2);
1034  uvc[0] = (info[0][0] + info[1][0] + info[2][0]) / 3.0;
1035  uvc[1] = (info[0][1] + info[1][1] + info[2][1]) / 3.0;
1036 
1037  Array<OneD, NekDouble> locc = m_cadsurf->P(uvc);
1038  NekDouble d4 = m_mesh->m_octree->Query(locc);
1039 
1040  NekDouble d = (d1 + d2 + d3 + d4) / 4.0;
1041 
1042  vector<bool> valid(3);
1043  valid[0] = r[0] < d * 1.5;
1044  valid[1] = r[1] < d * 1.5;
1045  valid[2] = r[2] < d * 1.5;
1046 
1047  vector<bool> angValid(3);
1048  angValid[0] = a[0] / M_PI * 180.0 > 20.0 && a[0] / M_PI * 180.0 < 120.0;
1049  angValid[1] = a[1] / M_PI * 180.0 > 20.0 && a[1] / M_PI * 180.0 < 120.0;
1050  angValid[2] = a[2] / M_PI * 180.0 > 20.0 && a[2] / M_PI * 180.0 < 120.0;
1051 
1052  int numValid = 0;
1053  int numAngValid = 0;
1054  for (int j = 0; j < 3; j++)
1055  {
1056  if (valid[j])
1057  {
1058  numValid++;
1059  }
1060  if (angValid[j])
1061  {
1062  numAngValid++;
1063  }
1064  }
1065 
1066  // if numvalid is zero no work to be done
1067  /*if (numValid != 3)
1068  {
1069  AddNewPoint(uvc);
1070  }*/
1071 
1072  if (numValid != 3 || numAngValid != 3)
1073  {
1074  // break the bad edge
1075  /*int a=0, b=0;
1076  if(!valid[0])
1077  {
1078  a = 0;
1079  b = 1;
1080  }
1081  else if(!valid[1])
1082  {
1083  a = 1;
1084  b = 2;
1085  }
1086  else if(!valid[2])
1087  {
1088  a = 2;
1089  b = 0;
1090  }
1091 
1092  Array<OneD, NekDouble> uvn(2);
1093  uvn[0] = (info[a][0] + info[b][0]) / 2.0;
1094  uvn[1] = (info[a][1] + info[b][1]) / 2.0;*/
1095  AddNewPoint(uvc);
1096  }
1097  }
1098 
1099  if (m_stienerpoints.size() == pointBefore)
1100  {
1101  return false;
1102  }
1103  else
1104  {
1105  return true;
1106  }
1107 }
1108 
1109 void FaceMesh::AddNewPoint(Array<OneD, NekDouble> uv)
1110 {
1111  // adds a new point but checks that there are no other points nearby first
1112  Array<OneD, NekDouble> np = m_cadsurf->P(uv);
1113  NekDouble npDelta = m_mesh->m_octree->Query(np);
1114 
1115  NodeSharedPtr n = boost::shared_ptr<Node>(
1116  new Node(m_mesh->m_numNodes++, np[0], np[1], np[2]));
1117 
1118  bool add = true;
1119 
1120  for (int i = 0; i < orderedLoops.size(); i++)
1121  {
1122  for (int j = 0; j < orderedLoops[i].size(); j++)
1123  {
1124  NekDouble r = orderedLoops[i][j]->Distance(n);
1125 
1126  if (r < npDelta / 2.0)
1127  {
1128  add = false;
1129  break;
1130  }
1131  }
1132  }
1133 
1134  if (add)
1135  {
1136  for (int i = 0; i < m_stienerpoints.size(); i++)
1137  {
1138  NekDouble r = m_stienerpoints[i]->Distance(n);
1139 
1140  if (r < npDelta / 2.0)
1141  {
1142  add = false;
1143  break;
1144  }
1145  }
1146  }
1147 
1148  if (add)
1149  {
1150  n->SetCADSurf(m_id, m_cadsurf, uv);
1151  m_stienerpoints.push_back(n);
1152  }
1153 }
1154 
1155 void FaceMesh::OrientateCurves()
1156 {
1157  // this could be a second run on orentate so clear some info
1158  orderedLoops.clear();
1159 
1160  // create list of bounding loop nodes
1161  for (int i = 0; i < m_edgeloops.size(); i++)
1162  {
1163  vector<NodeSharedPtr> cE;
1164  for (int j = 0; j < m_edgeloops[i]->edges.size(); j++)
1165  {
1166  int cid = m_edgeloops[i]->edges[j]->GetId();
1167  vector<NodeSharedPtr> edgePoints =
1168  m_curvemeshes[cid]->GetMeshPoints();
1169 
1170  int numPoints = m_curvemeshes[cid]->GetNumPoints();
1171 
1172  if (m_edgeloops[i]->edgeo[j] == CADOrientation::eForwards)
1173  {
1174  for (int k = 0; k < numPoints - 1; k++)
1175  {
1176  cE.push_back(edgePoints[k]);
1177  }
1178  }
1179  else
1180  {
1181  for (int k = numPoints - 1; k > 0; k--)
1182  {
1183  cE.push_back(edgePoints[k]);
1184  }
1185  }
1186  }
1187  orderedLoops.push_back(cE);
1188  }
1189 }
1190 }
1191 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Basic information about an element.
Definition: ElementConfig.h:50
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
Represents an edge which joins two points.
Definition: Edge.h:58
boost::shared_ptr< TriangleInterface > TriangleInterfaceSharedPtr
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
STL namespace.
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
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
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
boost::unordered_set< EdgeSharedPtr, EdgeHash > EdgeSet
Definition: Edge.h:162