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