Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
BLMesh.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: TetMesh.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: tet meshing methods
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <boost/geometry.hpp>
37 #include <boost/geometry/geometries/box.hpp>
38 #include <boost/geometry/geometries/point.hpp>
39 #include <boost/geometry/index/rtree.hpp>
40 
42 
45 
48 
49 #include <algorithm>
50 
51 namespace bg = boost::geometry;
52 namespace bgi = boost::geometry::index;
53 
54 typedef bg::model::point<double, 3, bg::cs::cartesian> point;
55 typedef bg::model::box<point> box;
56 typedef std::pair<box, unsigned int> boxI;
57 
58 using namespace std;
59 namespace Nektar
60 {
61 namespace NekMeshUtils
62 {
63 
65 {
66  NekDouble xmin = numeric_limits<double>::max(),
67  xmax = -1.0 * numeric_limits<double>::max(),
68  ymin = numeric_limits<double>::max(),
69  ymax = -1.0 * numeric_limits<double>::max(),
70  zmin = numeric_limits<double>::max(),
71  zmax = -1.0 * numeric_limits<double>::max();
72 
73  vector<NodeSharedPtr> ns = el->GetVertexList();
74  for (int i = 0; i < ns.size(); i++)
75  {
76  xmin = min(xmin, ns[i]->m_x);
77  xmax = max(xmax, ns[i]->m_x);
78  ymin = min(ymin, ns[i]->m_y);
79  ymax = max(ymax, ns[i]->m_y);
80  zmin = min(zmin, ns[i]->m_z);
81  zmax = max(zmax, ns[i]->m_z);
82  }
83 
84  return box(point(xmin - ov, ymin - ov, zmin - ov),
85  point(xmax + ov, ymax + ov, zmax + ov));
86 }
87 
88 inline box GetBox(vector<ElementSharedPtr> els, NekDouble ov)
89 {
90  NekDouble xmin = numeric_limits<double>::max(),
91  xmax = -1.0 * numeric_limits<double>::max(),
92  ymin = numeric_limits<double>::max(),
93  ymax = -1.0 * numeric_limits<double>::max(),
94  zmin = numeric_limits<double>::max(),
95  zmax = -1.0 * numeric_limits<double>::max();
96 
97  for (int j = 0; j < els.size(); j++)
98  {
99  vector<NodeSharedPtr> ns = els[j]->GetVertexList();
100  for (int i = 0; i < ns.size(); i++)
101  {
102  xmin = min(xmin, ns[i]->m_x);
103  xmax = max(xmax, ns[i]->m_x);
104  ymin = min(ymin, ns[i]->m_y);
105  ymax = max(ymax, ns[i]->m_y);
106  zmin = min(zmin, ns[i]->m_z);
107  zmax = max(zmax, ns[i]->m_z);
108  }
109  }
110 
111  return box(point(xmin - ov, ymin - ov, zmin - ov),
112  point(xmax + ov, ymax + ov, zmax + ov));
113 }
114 
116 {
117  return box(point(n->m_x - ov, n->m_y - ov, n->m_z - ov),
118  point(n->m_x + ov, n->m_y + ov, n->m_z + ov));
119 }
120 
122 {
123  Setup();
124 
125  BuildElements();
126 
127  GrowLayers();
128 
129  Shrink();
130 
132  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
133  {
134  vector<blInfoSharedPtr> infos = m_nToNInfo[bit->first];
135  for (int i = 0; i < infos.size(); i++)
136  {
137  if (bit->second->bl > infos[i]->bl + 1)
138  {
139  cout << "non smooth error " << bit->second->bl << " "
140  << infos[i]->bl << endl;
141  }
142  }
143  }
144 
145  for (int i = 0; i < m_mesh->m_element[3].size(); i++)
146  {
147  ElementSharedPtr el = m_mesh->m_element[3][i];
148  if (!IsPrismValid(el))
149  {
150  cout << "validity error " << el->GetId() << endl;
151  }
152  }
153 
154  /*m_mesh->m_element[2] = m_psuedoSurface;
155  m_mesh->m_element[3].clear();
156  m_mesh->m_expDim = 2;*/
157 }
158 
159 map<NodeSharedPtr, NodeSharedPtr> BLMesh::GetSymNodes()
160 {
161  map<NodeSharedPtr, NodeSharedPtr> ret;
162 
164  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
165  {
166  if (!bit->second->onSym)
167  {
168  continue;
169  }
170  CADSurfSharedPtr s = m_mesh->m_cad->GetSurf(bit->second->symsurf);
171  Array<OneD, NekDouble> loc = bit->second->pNode->GetLoc();
173  uv = s->locuv(loc);
174  bit->second->pNode->SetCADSurf(bit->second->symsurf, s, uv);
175  ret[bit->first] = bit->second->pNode;
176  }
177  return ret;
178 }
179 
181 {
182  vector<NodeSharedPtr> ns1 = el->GetVertexList();
183  Array<OneD, NekDouble> N1 = el->Normal(true);
184 
186  V[0] = n->m_x - ns1[0]->m_x;
187  V[1] = n->m_y - ns1[0]->m_y;
188  V[2] = n->m_z - ns1[0]->m_z;
189 
190  NekDouble Vmag = sqrt(V[0] * V[0] + V[1] * V[1] + V[2] * V[2]);
191 
192  NekDouble ang = (N1[0] * V[0] + N1[1] * V[1] + N1[2] * V[2]) / Vmag;
193 
194  return ang > 0.17;
195 }
196 
197 void BLMesh::GrowLayers()
198 {
200 
201  // setup up a tree which is formed of boxes of each surface plus some
202  // extra room (ideal bl thick)
203 
204  // in each iteration a tree is made for all the triangles in each surface
205  // when considering to stop a bounary layer growing, it first
206  // looks at the top tree to find surfaces which are canditates
207  // it then searches the subtrees for each triangle canditate
208  // it then does a distance calcation.
209  // if a boundary layer should be close to that from another surface, it
210  // should stop
211 
212  map<int, vector<ElementSharedPtr> > psElements;
213  for (int i = 0; i < m_mesh->m_element[2].size(); i++)
214  {
215  ElementSharedPtr el = m_mesh->m_element[2][i];
217  find(m_blsurfs.begin(), m_blsurfs.end(), el->m_parentCAD->GetId());
218 
220  m_symSurfs.begin(), m_symSurfs.end(), el->m_parentCAD->GetId());
221 
222  if (f == m_blsurfs.end() && s == m_symSurfs.end())
223  {
224  psElements[el->m_parentCAD->GetId()].push_back(el);
225  }
226  }
227  for (int i = 0; i < m_psuedoSurface.size(); i++)
228  {
229  psElements[m_psuedoSurface[i]->m_parentCAD->GetId()].push_back(
230  m_psuedoSurface[i]);
231  }
232 
233  bgi::rtree<boxI, bgi::quadratic<16> > TopTree;
234  map<int, bgi::rtree<boxI, bgi::quadratic<16> > > SubTrees;
235 
236  // ofstream file;
237  // file.open("pts.3D");
238  // file << "x y z value" << endl;
239 
240  for (int l = 1; l < m_layer; l++)
241  {
242  NekDouble delta = (m_layerT[l] - m_layerT[l - 1]);
243  TopTree.clear();
244  SubTrees.clear();
245  map<int, vector<ElementSharedPtr> >::iterator it;
246  for (it = psElements.begin(); it != psElements.end(); it++)
247  {
248  TopTree.insert(make_pair(GetBox(it->second, m_bl), it->first));
249  vector<boxI> toInsert;
250  for (int i = 0; i < it->second.size(); i++)
251  {
252  toInsert.push_back(make_pair(GetBox(it->second[i], m_bl), i));
253  }
254  SubTrees[it->first].insert(toInsert.begin(), toInsert.end());
255  }
256 
257  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
258  {
259  if (bit->second->stopped)
260  {
261  continue;
262  }
263 
264  vector<boxI> results;
265  TopTree.query(bgi::intersects(point(bit->second->pNode->m_x,
266  bit->second->pNode->m_y,
267  bit->second->pNode->m_z)),
268  back_inserter(results));
269  set<int> surfs;
270  for (int i = 0; i < results.size(); i++)
271  {
273  bit->second->surfs.find(results[i].second);
274  if (f == bit->second->surfs.end())
275  {
276  // hit
277  surfs.insert(results[i].second);
278  }
279  }
280 
281  set<int>::iterator iit;
282  bool hit = false;
283  for (iit = surfs.begin(); iit != surfs.end(); iit++)
284  {
285  results.clear();
286  SubTrees[*iit].query(
287  bgi::intersects(GetBox(bit->second->pNode, m_bl)),
288  back_inserter(results));
289  for (int i = 0; i < results.size(); i++)
290  {
291  if (Infont(bit->second->pNode,
292  psElements[*iit][results[i].second]))
293  {
294  NekDouble prox =
295  Proximity(bit->second->pNode,
296  psElements[*iit][results[i].second]);
297  if (prox < delta * 2.5)
298  {
299  hit = true;
300  // cout << "hit" << endl;
301  bit->second->stopped = true;
302  /*file << bit->first->m_x << " " << bit->first->m_y
303  << " " << bit->first->m_z << " " << l << endl;
304  file << bit->second->pNode->m_x << " " <<
305  bit->second->pNode->m_y << " " <<
306  bit->second->pNode->m_z << " " << l << endl;
307  m_mesh->m_element[2].clear();
308  m_mesh->m_expDim--;
309  m_mesh->m_element[2].push_back(psElements[*iit][results[i].second]);*/
310  break;
311  // return;
312  }
313  }
314  }
315  if (hit)
316  break;
317  }
318  }
319 
320  // if after proximity scanning all is okay, advance the layer
321  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
322  {
323  if (bit->second->stopped)
324  {
325  continue;
326  }
327 
328  // test the smoothness
329  bool shouldStop = false;
330  vector<blInfoSharedPtr> ne = m_nToNInfo[bit->first];
331  for (int i = 0; i < ne.size(); i++)
332  {
333  if (ne[i]->bl < bit->second->bl)
334  {
335  shouldStop = true;
336  break;
337  }
338  }
339  if (shouldStop)
340  {
341  bit->second->stopped = true;
342  continue;
343  }
344 
345  bit->second->AlignNode(m_layerT[l]);
346  bit->second->bl = l;
347  }
348  }
349  // file.close();
350 }
351 
352 inline bool sign(NekDouble a, NekDouble b)
353 {
354  return (a * b > 0.0);
355 }
356 
358 {
359  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
360 }
361 
362 NekDouble BLMesh::Proximity(NodeSharedPtr n, ElementSharedPtr el)
363 {
364  vector<NodeSharedPtr> ns = el->GetVertexList();
365  Array<OneD, NekDouble> B = ns[0]->GetLoc();
367  E0[0] = ns[1]->m_x - ns[0]->m_x;
368  E0[1] = ns[1]->m_y - ns[0]->m_y;
369  E0[2] = ns[1]->m_z - ns[0]->m_z;
371  E1[0] = ns[2]->m_x - ns[0]->m_x;
372  E1[1] = ns[2]->m_y - ns[0]->m_y;
373  E1[2] = ns[2]->m_z - ns[0]->m_z;
374  Array<OneD, NekDouble> P = n->GetLoc();
375 
376  NekDouble a = Dot(E0, E0);
377  NekDouble b = Dot(E0, E1);
378  NekDouble c = Dot(E1, E1);
379 
381  Vmath::Vsub(3, B, 1, P, 1, BP, 1);
382 
383  NekDouble d = Dot(E0, BP);
384  NekDouble e = Dot(E1, BP);
385 
386  NekDouble det = a * c - b * b;
387  NekDouble s = b * e - c * d, t = b * d - a * e;
388 
389  if (s + t <= det)
390  {
391  if (s < 0)
392  {
393  if (t < 0)
394  {
395  t = 0;
396  s = 0;
397  }
398  else
399  {
400  s = 0;
401  t = (e >= 0 ? 0 : (-e >= c ? 1 : -e / c));
402  }
403  }
404  else if (t < 0)
405  {
406  t = 0;
407  s = (d >= 0 ? 0 : (-d >= a ? 1 : -d / a));
408  }
409  else
410  {
411  NekDouble invDet = 1.0 / det;
412  s *= invDet;
413  t *= invDet;
414  }
415  }
416  else
417  {
418  if (s < 0)
419  {
420  s = 0;
421  t = 1;
422  }
423  else if (t < 0)
424  {
425  s = 1;
426  t = 0;
427  }
428  else
429  {
430  NekDouble numer = c + e - b - d;
431  if (numer <= 0)
432  {
433  s = 0;
434  }
435  else
436  {
437  NekDouble denom = a - 2 * b + c;
438  s = (numer >= denom ? 1 : numer / denom);
439  }
440  t = 1 - s;
441  }
442  }
443 
444  NodeSharedPtr point = boost::shared_ptr<Node>(
445  new Node(0, B[0] + s * E0[0] + t * E1[0], B[1] + s * E0[1] + t * E1[1],
446  B[2] + s * E0[2] + t * E1[2]));
447 
448  return n->Distance(point);
449 }
450 
451 bool BLMesh::TestIntersectionEl(ElementSharedPtr e1, ElementSharedPtr e2)
452 {
453  vector<NodeSharedPtr> ns1 = e1->GetVertexList();
454  vector<NodeSharedPtr> ns2 = e2->GetVertexList();
455  if (ns1[0] == ns2[0] || ns1[0] == ns2[1] || ns1[0] == ns2[2] ||
456  ns1[1] == ns2[0] || ns1[1] == ns2[1] || ns1[1] == ns2[2] ||
457  ns1[2] == ns2[0] || ns1[2] == ns2[1] || ns1[2] == ns2[2])
458  {
459  return false;
460  }
461 
462  Array<OneD, NekDouble> N1(3), N2(3);
463  NekDouble d1, d2;
464 
465  N1[0] = (ns1[1]->m_y - ns1[0]->m_y) * (ns1[2]->m_z - ns1[0]->m_z) -
466  (ns1[2]->m_y - ns1[0]->m_y) * (ns1[1]->m_z - ns1[0]->m_z);
467  N1[1] = -1.0 * ((ns1[1]->m_x - ns1[0]->m_x) * (ns1[2]->m_z - ns1[0]->m_z) -
468  (ns1[2]->m_x - ns1[0]->m_x) * (ns1[1]->m_z - ns1[0]->m_z));
469  N1[2] = (ns1[1]->m_x - ns1[0]->m_x) * (ns1[2]->m_y - ns1[0]->m_y) -
470  (ns1[2]->m_x - ns1[0]->m_x) * (ns1[1]->m_y - ns1[0]->m_y);
471 
472  N2[0] = (ns2[1]->m_y - ns2[0]->m_y) * (ns2[2]->m_z - ns2[0]->m_z) -
473  (ns2[2]->m_y - ns2[0]->m_y) * (ns2[1]->m_z - ns2[0]->m_z);
474  N2[1] = -1.0 * ((ns2[1]->m_x - ns2[0]->m_x) * (ns2[2]->m_z - ns2[0]->m_z) -
475  (ns2[2]->m_x - ns2[0]->m_x) * (ns2[1]->m_z - ns2[0]->m_z));
476  N2[2] = (ns2[1]->m_x - ns2[0]->m_x) * (ns2[2]->m_y - ns2[0]->m_y) -
477  (ns2[2]->m_x - ns2[0]->m_x) * (ns2[1]->m_y - ns2[0]->m_y);
478 
479  d1 = -1.0 *
480  (N1[0] * ns1[0]->m_x + N1[1] * ns1[0]->m_y + N1[2] * ns1[0]->m_z);
481  d2 = -1.0 *
482  (N2[0] * ns2[0]->m_x + N2[1] * ns2[0]->m_y + N2[2] * ns2[0]->m_z);
483 
484  Array<OneD, NekDouble> dv1(3), dv2(3);
485 
486  dv1[0] =
487  N2[0] * ns1[0]->m_x + N2[1] * ns1[0]->m_y + N2[2] * ns1[0]->m_z + d2;
488  dv1[1] =
489  N2[0] * ns1[1]->m_x + N2[1] * ns1[1]->m_y + N2[2] * ns1[1]->m_z + d2;
490  dv1[2] =
491  N2[0] * ns1[2]->m_x + N2[1] * ns1[2]->m_y + N2[2] * ns1[2]->m_z + d2;
492 
493  dv2[0] =
494  N1[0] * ns2[0]->m_x + N1[1] * ns2[0]->m_y + N1[2] * ns2[0]->m_z + d1;
495  dv2[1] =
496  N1[0] * ns2[1]->m_x + N1[1] * ns2[1]->m_y + N1[2] * ns2[1]->m_z + d1;
497  dv2[2] =
498  N1[0] * ns2[2]->m_x + N1[1] * ns2[2]->m_y + N1[2] * ns2[2]->m_z + d1;
499 
500  if (sign(dv1[0], dv1[1]) && sign(dv1[1], dv1[2]))
501  {
502  return false;
503  }
504  if (sign(dv2[0], dv2[1]) && sign(dv2[1], dv2[2]))
505  {
506  return false;
507  }
508 
510  D[0] = N1[1] * N2[2] - N1[2] * N2[1];
511  D[1] = -1.0 * (N1[0] * N2[2] - N1[2] * N2[0]);
512  D[2] = N1[0] * N2[1] - N1[0] * N2[1];
513 
514  int base1 = 0, base2 = 0;
515  if (!sign(dv2[0], dv2[1]) && sign(dv2[1], dv2[2]))
516  {
517  base2 = 0;
518  }
519  else if (!sign(dv2[1], dv2[2]) && sign(dv2[2], dv2[0]))
520  {
521  base2 = 1;
522  }
523  else if (!sign(dv2[2], dv2[0]) && sign(dv2[0], dv2[1]))
524  {
525  base2 = 2;
526  }
527  else
528  {
529  cout << "base not set" << endl;
530  }
531 
532  if (!sign(dv1[0], dv1[1]) && sign(dv1[1], dv1[2]))
533  {
534  base1 = 0;
535  }
536  else if (!sign(dv1[1], dv1[2]) && sign(dv1[2], dv1[0]))
537  {
538  base1 = 1;
539  }
540  else if (!sign(dv1[2], dv1[0]) && sign(dv1[0], dv1[1]))
541  {
542  base1 = 2;
543  }
544  else
545  {
546  cout << "base not set" << endl;
547  }
548 
549  Array<OneD, NekDouble> p1(3), p2(3);
550 
551  p1[0] = D[0] * ns1[0]->m_x + D[1] * ns1[0]->m_y + D[2] * ns1[0]->m_z;
552  p1[1] = D[0] * ns1[1]->m_x + D[1] * ns1[1]->m_y + D[2] * ns1[1]->m_z;
553  p1[2] = D[0] * ns1[2]->m_x + D[1] * ns1[2]->m_y + D[2] * ns1[2]->m_z;
554 
555  p2[0] = D[0] * ns2[0]->m_x + D[1] * ns2[0]->m_y + D[2] * ns2[0]->m_z;
556  p2[1] = D[0] * ns2[1]->m_x + D[1] * ns2[1]->m_y + D[2] * ns2[1]->m_z;
557  p2[2] = D[0] * ns2[2]->m_x + D[1] * ns2[2]->m_y + D[2] * ns2[2]->m_z;
558 
559  NekDouble t11, t12, t21, t22;
560  int o1 = 0, o2 = 0;
561  if (base1 == 0)
562  {
563  o1 = 1;
564  o2 = 2;
565  }
566  else if (base1 == 1)
567  {
568  o1 = 2;
569  o2 = 0;
570  }
571  else if (base1 == 2)
572  {
573  o1 = 0;
574  o2 = 1;
575  }
576 
577  t11 = p1[o1] + (p1[base1] - p1[o1]) * dv1[o1] / (dv1[o1] - dv1[base1]);
578  t12 = p1[o2] + (p1[base1] - p1[o2]) * dv1[o2] / (dv1[o2] - dv1[base1]);
579 
580  if (base2 == 0)
581  {
582  o1 = 1;
583  o2 = 2;
584  }
585  else if (base2 == 1)
586  {
587  o1 = 2;
588  o2 = 0;
589  }
590  else if (base2 == 2)
591  {
592  o1 = 0;
593  o2 = 1;
594  }
595 
596  t21 = p2[o1] + (p2[base2] - p2[o1]) * dv2[o1] / (dv2[o1] - dv2[base2]);
597  t22 = p2[o2] + (p2[base2] - p2[o2]) * dv2[o2] / (dv2[o2] - dv2[base2]);
598 
599  if (t11 > t12)
600  {
601  swap(t11, t12);
602  }
603  if (t21 > t22)
604  {
605  swap(t21, t22);
606  }
607 
608  if (t21 < t11)
609  {
610  swap(t11, t21);
611  swap(t12, t22);
612  }
613 
614  if (!sign(t21 - t11, t22 - t11) || !sign(t21 - t12, t22 - t12))
615  {
616  return true;
617  }
618 
619  return false;
620 }
621 
622 void BLMesh::Shrink()
623 {
625  bool smsh = true;
626  while (smsh)
627  {
628  smsh = false;
629 
630  vector<ElementSharedPtr> inv;
631  for (int i = 0; i < m_mesh->m_element[3].size(); i++)
632  {
633  ElementSharedPtr el = m_mesh->m_element[3][i];
634  if (!IsPrismValid(el))
635  {
636  inv.push_back(el);
637  }
638  }
639 
640  smsh = (inv.size() > 0);
641 
642  for (int i = 0; i < inv.size(); i++)
643  {
644  ElementSharedPtr t = m_priToTri[inv[i]];
645  vector<blInfoSharedPtr> bls;
646  vector<NodeSharedPtr> ns = t->GetVertexList();
647  for (int j = 0; j < ns.size(); j++)
648  {
649  bls.push_back(m_blData[ns[j]]);
650  }
651  bool repeat = true;
652  while (repeat)
653  {
654  repeat = false;
655  int mx = 0;
656  for (int j = 0; j < 3; j++)
657  {
658  mx = max(mx, bls[j]->bl);
659  }
660  ASSERTL0(mx > 0, "shrinking to nothing");
661  for (int j = 0; j < 3; j++)
662  {
663  if (bls[j]->bl < mx)
664  {
665  continue;
666  }
667  bls[j]->bl--;
668  bls[j]->AlignNode(m_layerT[bls[j]->bl]);
669  }
670  if (!IsPrismValid(inv[i]))
671  {
672  repeat = true;
673  }
674  }
675  }
676 
677  bool repeat = true;
678  while (repeat)
679  {
680  repeat = false;
681  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
682  {
683  vector<blInfoSharedPtr> infos = m_nToNInfo[bit->first];
684  for (int i = 0; i < infos.size(); i++)
685  {
686  if (bit->second->bl > infos[i]->bl + 1)
687  {
688  bit->second->bl--;
689  bit->second->AlignNode(m_layerT[bit->second->bl]);
690  repeat = true;
691  }
692  }
693  }
694  }
695  }
696 }
697 
698 bool BLMesh::IsPrismValid(ElementSharedPtr el)
699 {
700  NekDouble mn = numeric_limits<double>::max();
701  NekDouble mx = -1.0 * numeric_limits<double>::max();
702  vector<NodeSharedPtr> ns = el->GetVertexList();
703  NekVector<NekDouble> X(6), Y(6), Z(6);
704  for (int j = 0; j < ns.size(); j++)
705  {
706  X(j) = ns[j]->m_x;
707  Y(j) = ns[j]->m_y;
708  Z(j) = ns[j]->m_z;
709  }
710  NekVector<NekDouble> x1(6), y1(6), z1(6), x2(6), y2(6), z2(6), x3(6), y3(6),
711  z3(6);
712 
713  x1 = m_deriv[0] * X;
714  y1 = m_deriv[0] * Y;
715  z1 = m_deriv[0] * Z;
716  x2 = m_deriv[1] * X;
717  y2 = m_deriv[1] * Y;
718  z2 = m_deriv[1] * Z;
719  x3 = m_deriv[2] * X;
720  y3 = m_deriv[2] * Y;
721  z3 = m_deriv[2] * Z;
722 
723  for (int j = 0; j < 6; j++)
724  {
725  DNekMat dxdz(3, 3, 1.0, eFULL);
726  dxdz(0, 0) = x1(j);
727  dxdz(0, 1) = x2(j);
728  dxdz(0, 2) = x3(j);
729  dxdz(1, 0) = y1(j);
730  dxdz(1, 1) = y2(j);
731  dxdz(1, 2) = y3(j);
732  dxdz(2, 0) = z1(j);
733  dxdz(2, 1) = z2(j);
734  dxdz(2, 2) = z3(j);
735 
736  NekDouble jacDet =
737  dxdz(0, 0) * (dxdz(1, 1) * dxdz(2, 2) - dxdz(2, 1) * dxdz(1, 2)) -
738  dxdz(0, 1) * (dxdz(1, 0) * dxdz(2, 2) - dxdz(2, 0) * dxdz(1, 2)) +
739  dxdz(0, 2) * (dxdz(1, 0) * dxdz(2, 1) - dxdz(2, 0) * dxdz(1, 1));
740  mn = min(mn, jacDet);
741  mx = max(mx, jacDet);
742  }
743 
744  /*SpatialDomains::GeometrySharedPtr geom = el->GetGeom(3);
745  SpatialDomains::GeomFactorsSharedPtr gfac = geom->GetGeomFactors();
746 
747  cout << mn << " " << mx << " " << (mn > 0) << " " << gfac->IsValid() <<
748  endl;*/
749 
750  return mn > 0;
751 }
752 
753 void BLMesh::BuildElements()
754 {
755  // make prisms
756  map<CADOrientation::Orientation, vector<int> > baseTri;
757  map<CADOrientation::Orientation, vector<int> > topTri;
758 
759  vector<int> tmp;
760  // back-base
761  tmp.push_back(0);
762  tmp.push_back(4);
763  tmp.push_back(1);
764  baseTri[CADOrientation::eBackwards] = tmp;
765  tmp.clear();
766  // for-base
767  tmp.push_back(0);
768  tmp.push_back(1);
769  tmp.push_back(4);
770  baseTri[CADOrientation::eForwards] = tmp;
771  // back-top
772  tmp.clear();
773  tmp.push_back(3);
774  tmp.push_back(5);
775  tmp.push_back(2);
776  topTri[CADOrientation::eBackwards] = tmp;
777  // for-top
778  tmp.clear();
779  tmp.push_back(3);
780  tmp.push_back(2);
781  tmp.push_back(5);
782  topTri[CADOrientation::eForwards] = tmp;
783 
784  ElmtConfig pconf(LibUtilities::ePrism, 1, false, false);
785  ElmtConfig tconf(LibUtilities::eTriangle, 1, false, false);
786 
787  for (int i = 0; i < m_mesh->m_element[2].size(); i++)
788  {
789  ElementSharedPtr el = m_mesh->m_element[2][i];
791  find(m_blsurfs.begin(), m_blsurfs.end(), el->m_parentCAD->GetId());
792 
793  if (f == m_blsurfs.end())
794  {
795  // if this triangle is not in bl surfs continue
796  continue;
797  }
798 
799  vector<NodeSharedPtr> tn(3); // nodes for pseduo surface
800  vector<NodeSharedPtr> pn(6); // all prism nodes
801  vector<NodeSharedPtr> n = el->GetVertexList();
802  CADOrientation::Orientation o = el->m_parentCAD->Orientation();
803 
804  for (int j = 0; j < 3; j++)
805  {
806  pn[baseTri[o][j]] = n[j];
807  pn[topTri[o][j]] = m_blData[n[j]]->pNode;
808  tn[j] = m_blData[n[j]]->pNode;
809  }
810 
811  vector<int> tags;
812  tags.push_back(m_id);
814  LibUtilities::ePrism, pconf, pn, tags);
815  E->SetId(i);
816 
817  m_mesh->m_element[3].push_back(E);
818 
819  // tag of this element doesnt matter so can just be 1
821  LibUtilities::eTriangle, tconf, tn, tags);
822  m_psuedoSurface.push_back(T);
823 
824  T->m_parentCAD = el->m_parentCAD;
825 
826  m_priToTri[E] = el;
827  }
828 }
829 
830 NekDouble BLMesh::Visability(vector<ElementSharedPtr> tris,
832 {
833  NekDouble mn = numeric_limits<double>::max();
834 
835  for (int i = 0; i < tris.size(); i++)
836  {
837  Array<OneD, NekDouble> tmp = tris[i]->Normal(true);
838  NekDouble dt = tmp[0] * N[0] + tmp[1] * N[1] + tmp[2] * N[2];
839  mn = min(mn, dt);
840  }
841  return mn;
842 }
843 
844 Array<OneD, NekDouble> BLMesh::GetNormal(vector<ElementSharedPtr> tris)
845 {
846  // compile list of normals
847  vector<Array<OneD, NekDouble> > N;
848  for (int i = 0; i < tris.size(); i++)
849  {
850  N.push_back(tris[i]->Normal(true));
851  }
852 
853  vector<NekDouble> w(N.size());
854  Array<OneD, NekDouble> Np(3, 0.0);
855 
856  for (int i = 0; i < N.size(); i++)
857  {
858  w[i] = 1.0 / N.size();
859  }
860 
861  for (int i = 0; i < N.size(); i++)
862  {
863  Np[0] += w[i] * N[i][0];
864  Np[1] += w[i] * N[i][1];
865  Np[2] += w[i] * N[i][2];
866  }
867  NekDouble mag = sqrt(Np[0] * Np[0] + Np[1] * Np[1] + Np[2] * Np[2]);
868  Np[0] /= mag;
869  Np[1] /= mag;
870  Np[2] /= mag;
871 
872  Array<OneD, NekDouble> Ninital = Np;
873 
874  NekDouble dot = 0.0;
875  int ct = 0;
876  vector<NekDouble> a(N.size());
877  while (fabs(dot - 1) > 1e-6)
878  {
879  ct++;
880  Array<OneD, NekDouble> Nplast(3);
881  Nplast[0] = Np[0];
882  Nplast[1] = Np[1];
883  Nplast[2] = Np[2];
884 
885  NekDouble aSum = 0.0;
886  for (int i = 0; i < N.size(); i++)
887  {
888  NekDouble dot2 =
889  Np[0] * N[i][0] + Np[1] * N[i][1] + Np[2] * N[i][2];
890  if (fabs(dot2 - 1) < 1e-9)
891  {
892  a[i] = dot2 / fabs(dot2) * 1e-9;
893  }
894  else
895  {
896  a[i] = acos(dot2);
897  }
898 
899  aSum += a[i];
900  }
901 
902  NekDouble wSum = 0.0;
903  for (int i = 0; i < N.size(); i++)
904  {
905  w[i] = w[i] * a[i] / aSum;
906  wSum += w[i];
907  }
908 
909  for (int i = 0; i < N.size(); i++)
910  {
911  w[i] = w[i] / wSum;
912  }
913 
914  Array<OneD, NekDouble> NpN(3, 0.0);
915  for (int i = 0; i < N.size(); i++)
916  {
917  NpN[0] += w[i] * N[i][0];
918  NpN[1] += w[i] * N[i][1];
919  NpN[2] += w[i] * N[i][2];
920  }
921  mag = sqrt(NpN[0] * NpN[0] + NpN[1] * NpN[1] + NpN[2] * NpN[2]);
922  NpN[0] /= mag;
923  NpN[1] /= mag;
924  NpN[2] /= mag;
925 
926  Np[0] = 0.8 * NpN[0] + (1.0 - 0.8) * Np[0];
927  Np[1] = 0.8 * NpN[1] + (1.0 - 0.8) * Np[1];
928  Np[2] = 0.8 * NpN[2] + (1.0 - 0.8) * Np[2];
929  mag = sqrt(Np[0] * Np[0] + Np[1] * Np[1] + Np[2] * Np[2]);
930  Np[0] /= mag;
931  Np[1] /= mag;
932  Np[2] /= mag;
933 
934  dot = Np[0] * Nplast[0] + Np[1] * Nplast[1] + Np[2] * Nplast[2];
935 
936  if (ct > 100000)
937  {
938  cout << "run out of iterations" << endl;
939  Np = Ninital;
940  break;
941  }
942  }
943 
944  return Np;
945 }
946 
947 void BLMesh::Setup()
948 {
949  NekDouble a = 2.0 * (1.0 - m_prog) / (1.0 - pow(m_prog, m_layer + 1));
950  m_layerT = Array<OneD, NekDouble>(m_layer);
951  m_layerT[0] = a * m_prog * m_bl;
952  for (int i = 1; i < m_layer; i++)
953  {
954  m_layerT[i] = m_layerT[i - 1] + a * pow(m_prog, i) * m_bl;
955  }
956 
957  if(m_mesh->m_verbose)
958  {
959  cout << "First layer height " << m_layerT[0] << endl;
960  }
961 
962  // this sets up all the boundary layer normals data holder
963  set<int> symSurfs;
965  int ct = 0;
966  int failed = 0;
967 
968  // ofstream file1;
969  // file1.open("pts.3D");
970  // file1 << "X Y Z value" << endl;
971  for (it = m_mesh->m_vertexSet.begin(); it != m_mesh->m_vertexSet.end();
972  it++, ct++)
973  {
974  vector<pair<int, CADSurfSharedPtr> > ss = (*it)->GetCADSurfs();
975  vector<unsigned int> surfs;
976  for (int i = 0; i < ss.size(); i++)
977  {
978  surfs.push_back(ss[i].first);
979  }
980  sort(surfs.begin(), surfs.end());
981  vector<unsigned int> inter, diff;
982 
983  set_intersection(m_blsurfs.begin(), m_blsurfs.end(), surfs.begin(),
984  surfs.end(), back_inserter(inter));
985  set_symmetric_difference(inter.begin(), inter.end(), surfs.begin(),
986  surfs.end(), back_inserter(diff));
987 
988  // is somewhere on a bl surface
989  if (inter.size() > 0)
990  {
991  // initialise a new bl boudnary node
992  blInfoSharedPtr bln = boost::shared_ptr<blInfo>(new blInfo);
993  bln->oNode = (*it);
994  bln->stopped = false;
995 
996  // file1 << (*it)->m_x << " " << (*it)->m_y << " " << (*it)->m_z <<
997  // " " << ss.size() << endl;
998 
999  if (diff.size() > 0)
1000  {
1001  // if the diff size is greater than 1 there is a curve that
1002  // needs remeshing
1003  ASSERTL0(diff.size() <= 1, "not setup for curve bl refinement");
1004  symSurfs.insert(diff[0]);
1005  bln->symsurf = diff[0];
1006  bln->onSym = true;
1007  }
1008  else
1009  {
1010  bln->onSym = false;
1011  }
1012 
1013  m_blData[(*it)] = bln;
1014  }
1015  }
1016  // file1.close();
1017 
1018  // need a map from vertex idx to surface elements
1019  // but do not care about triangles which are not in the bl
1020  for (int i = 0; i < m_mesh->m_element[2].size(); i++)
1021  {
1023  find(m_blsurfs.begin(), m_blsurfs.end(),
1024  m_mesh->m_element[2][i]->m_parentCAD->GetId());
1025 
1026  if (f == m_blsurfs.end())
1027  {
1028  // if this triangle is not in bl surfs continue
1029  continue;
1030  }
1031 
1032  vector<NodeSharedPtr> ns = m_mesh->m_element[2][i]->GetVertexList();
1033  for (int j = 0; j < ns.size(); j++)
1034  {
1035  m_blData[ns[j]]->els.push_back(m_mesh->m_element[2][i]);
1036  m_blData[ns[j]]->surfs.insert(
1037  m_mesh->m_element[2][i]->m_parentCAD->GetId());
1038  }
1039  }
1040 
1042  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
1043  {
1044  // calculate mesh normal
1045  bit->second->N = GetNormal(bit->second->els);
1046 
1047  if (Visability(bit->second->els, bit->second->N) < 0.0)
1048  {
1049  cerr << "failed " << bit->first->m_x << " " << bit->first->m_y
1050  << " " << bit->first->m_z << " "
1051  << Visability(bit->second->els, bit->second->N) << endl;
1052  failed++;
1053  }
1054 
1055  Array<OneD, NekDouble> loc = bit->first->GetLoc();
1056  for (int k = 0; k < 3; k++)
1057  {
1058  loc[k] += bit->second->N[k] * m_layerT[0];
1059  }
1060 
1061  bit->second->pNode = boost::shared_ptr<Node>(
1062  new Node(m_mesh->m_numNodes++, loc[0], loc[1], loc[2]));
1063  bit->second->bl = 0;
1064  }
1065 
1066  m_symSurfs = vector<unsigned int>(symSurfs.begin(), symSurfs.end());
1067 
1068  // now need to enforce that all symmetry plane nodes have their normal
1069  // forced onto the symmetry surface
1070  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
1071  {
1072  if (!bit->second->onSym)
1073  {
1074  continue;
1075  }
1076 
1077  Array<OneD, NekDouble> uv(2);
1078  Array<OneD, NekDouble> loc = bit->second->pNode->GetLoc();
1079  m_mesh->m_cad->GetSurf(bit->second->symsurf)->ProjectTo(loc, uv);
1080 
1082  m_mesh->m_cad->GetSurf(bit->second->symsurf)->P(uv);
1083 
1085  N[0] = nl[0] - bit->first->m_x;
1086  N[1] = nl[1] - bit->first->m_y;
1087  N[2] = nl[2] - bit->first->m_z;
1088 
1089  NekDouble mag = sqrt(N[0] * N[0] + N[1] * N[1] + N[2] * N[2]);
1090  N[0] /= mag;
1091  N[1] /= mag;
1092  N[2] /= mag;
1093 
1094  bit->second->N = N;
1095  bit->second->AlignNode(m_layerT[0]);
1096  }
1097 
1098  // now smooth all the normals by distance weighted average
1099  // keep normals on curves constant
1100  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
1101  {
1102  set<int> added;
1103  added.insert(bit->first->m_id);
1104  for (int i = 0; i < bit->second->els.size(); i++)
1105  {
1106  vector<NodeSharedPtr> ns = bit->second->els[i]->GetVertexList();
1107  for (int j = 0; j < ns.size(); j++)
1108  {
1109  set<int>::iterator t = added.find(ns[j]->m_id);
1110  if (t == added.end())
1111  {
1112  m_nToNInfo[bit->first].push_back(m_blData[ns[j]]);
1113  }
1114  }
1115  }
1116  }
1117 
1118  for (int l = 0; l < 10; l++)
1119  {
1120  for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
1121  {
1122  if (bit->first->GetNumCADSurf() > 1)
1123  {
1124  continue;
1125  }
1126 
1127  Array<OneD, NekDouble> sumV(3, 0.0);
1128  vector<blInfoSharedPtr> data = m_nToNInfo[bit->first];
1129  NekDouble Dtotal = 0.0;
1130  for (int i = 0; i < data.size(); i++)
1131  {
1132  NekDouble d = bit->first->Distance(data[i]->oNode);
1133  Dtotal += d;
1134  sumV[0] += data[i]->N[0] / d;
1135  sumV[1] += data[i]->N[1] / d;
1136  sumV[2] += data[i]->N[2] / d;
1137  }
1138  sumV[0] *= Dtotal;
1139  sumV[1] *= Dtotal;
1140  sumV[2] *= Dtotal;
1141  NekDouble mag =
1142  sqrt(sumV[0] * sumV[0] + sumV[1] * sumV[1] + sumV[2] * sumV[2]);
1143  sumV[0] /= mag;
1144  sumV[1] /= mag;
1145  sumV[2] /= mag;
1146 
1148 
1149  N[0] = (1.0 - 0.8) * bit->second->N[0] + 0.8 * sumV[0];
1150  N[1] = (1.0 - 0.8) * bit->second->N[1] + 0.8 * sumV[1];
1151  N[2] = (1.0 - 0.8) * bit->second->N[2] + 0.8 * sumV[2];
1152 
1153  mag = sqrt(N[0] * N[0] + N[1] * N[1] + N[2] * N[2]);
1154  N[0] /= mag;
1155  N[1] /= mag;
1156  N[2] /= mag;
1157 
1158  bit->second->N = N;
1159  bit->second->AlignNode(m_layerT[0]);
1160  }
1161  }
1162 
1163  /*ofstream file;
1164  file.open("bl.lines");
1165  for(bit = m_blData.begin(); bit != m_blData.end(); bit++)
1166  {
1167  NekDouble l = 0.05;
1168  file << bit->first->m_x << ", " << bit->first->m_y << ", " <<
1169  bit->first->m_z << endl;
1170  file << bit->first->m_x + bit->second->N[0]*l << ", "
1171  << bit->first->m_y + bit->second->N[1]*l << ", "
1172  << bit->first->m_z + bit->second->N[2]*l << endl;
1173  file << endl;
1174  }
1175  file.close();*/
1176 
1177  ASSERTL0(failed == 0, "some normals failed to generate");
1178 
1180 
1181  Array<OneD, NekDouble> u1, v1, w1;
1182  LibUtilities::PointsManager()[pkey1]->GetPoints(u1, v1, w1);
1183 
1184  LibUtilities::NodalUtilPrism nodalPrism(1, u1, v1, w1);
1185 
1186  NekMatrix<NekDouble> Vandermonde = *nodalPrism.GetVandermonde();
1187  NekMatrix<NekDouble> VandermondeI = Vandermonde;
1188  VandermondeI.Invert();
1189 
1190  m_deriv[0] = *nodalPrism.GetVandermondeForDeriv(0) * VandermondeI;
1191  m_deriv[1] = *nodalPrism.GetVandermondeForDeriv(1) * VandermondeI;
1192  m_deriv[2] = *nodalPrism.GetVandermondeForDeriv(2) * VandermondeI;
1193 }
1194 }
1195 }
static NekDouble ang(NekDouble x, NekDouble y)
#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
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:27
STL namespace.
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
SharedMatrix GetVandermondeForDeriv(int dir)
Return the Vandermonde matrix of the derivative of the basis functions for the nodal distribution...
Definition: NodalUtil.cpp:135
bg::model::box< point > box
Definition: BLMesh.cpp:55
bool Infont(NodeSharedPtr n, ElementSharedPtr el)
Definition: BLMesh.cpp:180
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
PointsManagerT & PointsManager(void)
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
box GetBox(NodeSharedPtr n, NekDouble ov)
Definition: BLMesh.cpp:115
Specialisation of the NodalUtil class to support nodal prismatic elements.
Definition: NodalUtil.h:264
bg::model::point< double, 3, bg::cs::cartesian > point
Definition: BLMesh.cpp:54
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
Definition: Vmath.cpp:914
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:343
boost::shared_ptr< CADSurf > CADSurfSharedPtr
Definition: CADSurf.h:172
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::pair< box, unsigned int > boxI
Definition: BLMesh.cpp:56
SharedMatrix GetVandermonde()
Return the Vandermonde matrix for the nodal distribution.
Definition: NodalUtil.cpp:102
boost::shared_ptr< blInfo > blInfoSharedPtr
Definition: BLMesh.h:105
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:316
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
3D electrostatically spaced points on a Prism
Definition: PointsType.h:76