Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Octree.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: octree.h
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: cad object methods.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 //#include <algorithm>
37 //#include <limits>
38 
41 
42 using namespace std;
43 namespace Nektar
44 {
45 namespace NekMeshUtils
46 {
47 
49 {
50  // starting at master octant 0 move through succsesive m_octants which
51  // contain
52  // the point loc until a leaf is found
53  OctantSharedPtr n = m_masteroct;
54  int quad;
55 
56  bool found = false;
57 
58  while (!found)
59  {
60  Array<OneD, NekDouble> octloc = n->GetLoc();
61 
62  if (!(loc[0] < octloc[0]) && // forward
63  !(loc[1] < octloc[1]) && // forward
64  !(loc[2] < octloc[2])) // forward
65  {
66  quad = 0;
67  }
68  else if (!(loc[0] < octloc[0]) && // forward
69  !(loc[1] < octloc[1]) && // forward
70  !(loc[2] > octloc[2])) // back
71  {
72  quad = 1;
73  }
74  else if (!(loc[0] < octloc[0]) && // forward
75  !(loc[1] > octloc[1]) && // back
76  !(loc[2] < octloc[2])) // forward
77  {
78  quad = 2;
79  }
80  else if (!(loc[0] < octloc[0]) && // forward
81  !(loc[1] > octloc[1]) && // back
82  !(loc[2] > octloc[2])) // back
83  {
84  quad = 3;
85  }
86  else if (!(loc[0] > octloc[0]) && // back
87  !(loc[1] < octloc[1]) && // forward
88  !(loc[2] < octloc[2])) // forward
89  {
90  quad = 4;
91  }
92  else if (!(loc[0] > octloc[0]) && // back
93  !(loc[1] < octloc[1]) && // forward
94  !(loc[2] > octloc[2])) // back
95  {
96  quad = 5;
97  }
98  else if (!(loc[0] > octloc[0]) && // back
99  !(loc[1] > octloc[1]) && // back
100  !(loc[2] < octloc[2])) // forward
101  {
102  quad = 6;
103  }
104  else if (!(loc[0] > octloc[0]) && // back
105  !(loc[1] > octloc[1]) && // back
106  !(loc[2] > octloc[2])) // back
107  {
108  quad = 7;
109  }
110  else
111  {
112  ASSERTL0(false, "Cannot locate quadrant");
113  }
114 
115  n = n->GetChild(quad);
116 
117  if (n->IsLeaf())
118  {
119  found = true;
120  }
121  }
122  return n->GetDelta();
123 }
124 
125 void Octree::GetOctreeMesh(MeshSharedPtr m)
126 {
127  for (int i = 0; i < m_octants.size(); i++)
128  {
129  /*if(m_octants[i]->GetLocation() != eOnBoundary)
130  {
131  continue;
132  }*/
133 
134  vector<NodeSharedPtr> ns(8);
135 
136  ns[0] = boost::shared_ptr<Node>(new Node(0,
137  m_octants[i]->FX(eBack),
138  m_octants[i]->FX(eDown),
139  m_octants[i]->FX(eRight)));
140 
141  ns[1] = boost::shared_ptr<Node>(new Node(0,
142  m_octants[i]->FX(eForward),
143  m_octants[i]->FX(eDown),
144  m_octants[i]->FX(eRight)));
145 
146  ns[2] = boost::shared_ptr<Node>(new Node(0,
147  m_octants[i]->FX(eForward),
148  m_octants[i]->FX(eUp),
149  m_octants[i]->FX(eRight)));
150 
151  ns[3] = boost::shared_ptr<Node>(new Node(0,
152  m_octants[i]->FX(eBack),
153  m_octants[i]->FX(eUp),
154  m_octants[i]->FX(eRight)));
155 
156  ns[4] = boost::shared_ptr<Node>(new Node(0,
157  m_octants[i]->FX(eBack),
158  m_octants[i]->FX(eDown),
159  m_octants[i]->FX(eLeft)));
160 
161  ns[5] = boost::shared_ptr<Node>(new Node(0,
162  m_octants[i]->FX(eForward),
163  m_octants[i]->FX(eDown),
164  m_octants[i]->FX(eLeft)));
165 
166  ns[6] = boost::shared_ptr<Node>(new Node(0,
167  m_octants[i]->FX(eForward),
168  m_octants[i]->FX(eUp),
169  m_octants[i]->FX(eLeft)));
170 
171  ns[7] = boost::shared_ptr<Node>(new Node(0,
172  m_octants[i]->FX(eBack),
173  m_octants[i]->FX(eUp),
174  m_octants[i]->FX(eLeft)));
175 
176  vector<int> tags;
177  tags.push_back(0);
178  ElmtConfig conf(LibUtilities::eHexahedron, 1, false, false);
180  LibUtilities::eHexahedron, conf, ns, tags);
181  m->m_element[3].push_back(E);
182  }
183 }
184 
185 void Octree::Build()
186 {
187  Array<OneD, NekDouble> boundingBox = m_cad->GetBoundingBox();
188 
189  if (m_verbose)
190  cout << endl << "Octree system" << endl;
191 
192  // build curvature samples
193  CompileCuravturePointList();
194 
195  if (m_verbose)
196  cout << "\tCurvature samples: " << m_cpList.size() << endl;
197 
198  m_dim = max((boundingBox[1] - boundingBox[0]) / 2.0,
199  (boundingBox[3] - boundingBox[2]) / 2.0);
200 
201  m_dim = max(m_dim, (boundingBox[5] - boundingBox[4]) / 2.0);
202 
203  m_centroid = Array<OneD, NekDouble>(3);
204  m_centroid[0] = (boundingBox[1] + boundingBox[0]) / 2.0;
205  m_centroid[1] = (boundingBox[3] + boundingBox[2]) / 2.0;
206  m_centroid[2] = (boundingBox[5] + boundingBox[4]) / 2.0;
207 
208  // make master octant based on the bounding box of the domain
210  0, m_centroid[0], m_centroid[1], m_centroid[2], m_dim, m_cpList);
211 
212  SubDivide();
213 
214  m_octants.clear();
215  m_masteroct->CompileLeaves(m_octants);
216 
217  if (m_verbose)
218  {
219  cout << "\tOctants: " << m_octants.size() << endl;
220  }
221 
222  SmoothSurfaceOctants();
223 
224  PropagateDomain();
225 
226  SmoothAllOctants();
227 
228  if (m_verbose)
229  {
230  int elem = CountElemt();
231  cout << "\tPredicted mesh: " << elem << endl;
232  }
233 }
234 
235 void Octree::SubDivide()
236 {
237  bool repeat;
238  int ct = 1;
239 
240  m_numoct = 1;
241  m_masteroct->Subdivide(m_masteroct, m_numoct);
242 
243  if (m_verbose)
244  {
245  cout << "\tSubdivide iteration: ";
246  }
247 
248  do
249  {
250  ct++;
251  if (m_verbose)
252  {
253  cout << ct << " ";
254  cout.flush();
255  }
256  repeat = false;
257  m_octants.clear();
258  // grab a list of the leaves curently in the octree
259  m_masteroct->CompileLeaves(m_octants);
260 
261  VerifyNeigbours();
262 
263  // neeed to create a divide list, in the first list will be m_octants
264  // which need to
265  // subdivide based on curvature,
266  // in the next list will be ocants which need to subdivide to make sure
267  // level criteria is statisified for the previous list and so on
268  // the list will keep building till no more lists are required.
269  // the list will then be iterated through backwards to subdivide all the
270  // sublists in turn.
271  vector<vector<OctantSharedPtr> > dividelist;
272  set<int> inlist;
273  // build initial list
274  {
275  vector<OctantSharedPtr> sublist;
276  for (int i = 0; i < m_octants.size(); i++)
277  {
278  if (m_octants[i]->NeedDivide() &&
279  m_octants[i]->DX() / 4.0 > m_minDelta)
280  {
281  sublist.push_back(m_octants[i]);
282  inlist.insert(m_octants[i]->GetId());
283  repeat = true; // if there is a subdivision this whole
284  // process needs to be repeated
285  }
286  }
287  dividelist.push_back(sublist);
288  }
289  // then loop over building sublists until no more are required
290  int ct2 = 0;
291  while (true)
292  {
293  ct2++;
294  vector<OctantSharedPtr> newsublist,
295  previouslist = dividelist.back();
296  for (int i = 0; i < previouslist.size(); i++)
297  {
298  map<OctantFace, vector<OctantSharedPtr> > nlist =
299  previouslist[i]->GetNeigbours();
300  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
301  for (it = nlist.begin(); it != nlist.end(); it++)
302  {
303  for (int j = 0; j < it->second.size(); j++)
304  {
305  if (previouslist[i]->DX() < it->second[j]->DX())
306  {
308  inlist.find(it->second[j]->GetId());
309  if (s == inlist.end())
310  {
311  inlist.insert(it->second[j]->GetId());
312  newsublist.push_back(it->second[j]);
313  }
314  }
315  }
316  }
317  }
318  if (newsublist.size() == 0)
319  {
320  break;
321  }
322  else
323  {
324  dividelist.push_back(newsublist);
325  }
326  }
327 
328  vector<vector<OctantSharedPtr> >::reverse_iterator rit;
329  for (rit = dividelist.rbegin(); rit != dividelist.rend(); rit++)
330  {
331  vector<OctantSharedPtr> currentlist = *rit;
332  for (int i = 0; i < currentlist.size(); i++)
333  {
334  currentlist[i]->Subdivide(currentlist[i], m_numoct);
335  }
336  }
337  } while (repeat);
338 
339  if (m_verbose)
340  {
341  cout << endl;
342  }
343 }
344 
345 bool Octree::VerifyNeigbours()
346 {
347  // check all neibours
348  bool error = false;
349  for (int i = 0; i < m_octants.size(); i++)
350  {
351  bool valid = true;
352  map<OctantFace, vector<OctantSharedPtr> > nlist =
353  m_octants[i]->GetNeigbours();
354  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
355  for (it = nlist.begin(); it != nlist.end(); it++)
356  {
357  if (it->second.size() == 0)
358  {
359  NekDouble expectedfx;
360  switch (it->first)
361  {
362  case eUp:
363  expectedfx = m_centroid[1] + m_dim;
364  break;
365  case eDown:
366  expectedfx = m_centroid[1] - m_dim;
367  break;
368  case eLeft:
369  expectedfx = m_centroid[2] + m_dim;
370  break;
371  case eRight:
372  expectedfx = m_centroid[2] - m_dim;
373  break;
374  case eForward:
375  expectedfx = m_centroid[0] + m_dim;
376  break;
377  case eBack:
378  expectedfx = m_centroid[0] - m_dim;
379  break;
380  }
381  if (fabs(m_octants[i]->FX(it->first) - expectedfx) > 1E-6)
382  {
383  valid = false;
384  cout << "wall neigbour error" << endl;
385  cout << expectedfx << " " << m_octants[i]->FX(it->first)
386  << " " << it->first << endl;
387  }
388  }
389  else if (it->second.size() == 1)
390  {
391  if (!(m_octants[i]->DX() == it->second[0]->DX() ||
392  it->second[0]->DX() == 2.0 * m_octants[i]->DX()))
393  {
394  valid = false;
395  cout << " 1 neigbour error" << endl;
396  cout << m_octants[i]->DX() << " " << it->second[0]->DX()
397  << endl;
398  }
399  }
400  else if (it->second.size() == 4)
401  {
402  if (!(m_octants[i]->DX() / 2.0 == it->second[0]->DX()))
403  {
404  valid = false;
405  cout << "4 neibour error" << endl;
406  cout << m_octants[i]->DX() << " " << it->second[0]->DX()
407  << endl;
408  }
409  }
410  }
411  if (!valid)
412  {
413  error = true;
414  cout << "invalid neigbour config" << endl;
415  }
416  }
417  return !error;
418 }
419 
420 void Octree::SmoothSurfaceOctants()
421 {
422  // for all the m_octants which are surface containing and know their delta
423  // specification, look over all neighbours and ensure the specification
424  // between them is smooth
425  int ct = 0;
426 
427  do
428  {
429  ct = 0;
430  for (int i = 0; i < m_octants.size(); i++)
431  {
432  OctantSharedPtr oct = m_octants[i];
433 
434  if (oct->IsDeltaKnown())
435  {
436  vector<OctantSharedPtr> check;
437  map<OctantFace, vector<OctantSharedPtr> > nList =
438  oct->GetNeigbours();
439  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
440 
441  for (it = nList.begin(); it != nList.end(); it++)
442  {
443  for (int j = 0; j < it->second.size(); j++)
444  {
445  if (it->second[j]->IsDeltaKnown() &&
446  it->second[j]->GetDelta() < oct->GetDelta() &&
447  ddx(oct, it->second[j]) > 0.1)
448  {
449  check.push_back(it->second[j]);
450  }
451  }
452  }
453 
454  // for each neighbour listed in check_id, figure out the
455  // smoothed
456  // delta, and asign the miminum of these to nodes[i].GetDelta()
457  if (check.size() > 0)
458  {
459  NekDouble deltaSM = numeric_limits<double>::max();
460  for (int j = 0; j < check.size(); j++)
461  {
462  NekDouble r = oct->Distance(check[j]);
463 
464  if (0.099 * r + check[j]->GetDelta() < deltaSM)
465  {
466  deltaSM = 0.099 * r + check[j]->GetDelta();
467  }
468  }
469  oct->SetDelta(deltaSM);
470  ct += 1;
471  }
472  }
473  }
474  } while (ct > 0);
475 }
476 
477 void Octree::PropagateDomain()
478 {
479  // until all m_octants know their delta specifcation and orientaion
480  // look over all m_octants and if their neighours know either their
481  // orientation
482  // or specifcation calculate one for this octant
483  int ct = 0;
484 
485  do
486  {
487  ct = 0;
488  for (int i = 0; i < m_octants.size(); i++)
489  {
490  OctantSharedPtr oct = m_octants[i];
491 
492  if (!oct->IsDeltaKnown())
493  { // if delta has not been asigned
494  vector<OctantSharedPtr> known;
495  map<OctantFace, vector<OctantSharedPtr> > nList =
496  oct->GetNeigbours();
497  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
498 
499  for (it = nList.begin(); it != nList.end(); it++)
500  {
501  for (int j = 0; j < it->second.size(); j++)
502  {
503  if (it->second[j]->IsDeltaKnown())
504  {
505  known.push_back(it->second[j]);
506  }
507  }
508  }
509 
510  if (known.size() > 0)
511  {
512  vector<NekDouble> deltaPrime;
513  for (int j = 0; j < known.size(); j++)
514  {
515  NekDouble r = oct->Distance(known[j]);
516 
517  if (0.14 * r + known[j]->GetDelta() < m_maxDelta)
518  {
519  deltaPrime.push_back(0.14 * r +
520  known[j]->GetDelta());
521  }
522  else
523  {
524  deltaPrime.push_back(m_maxDelta);
525  }
526  }
527  NekDouble min = numeric_limits<double>::max();
528  for (int j = 0; j < deltaPrime.size(); j++)
529  {
530  if (deltaPrime[j] < min)
531  {
532  min = deltaPrime[j];
533  }
534  }
535  oct->SetDelta(min);
536  ct += 1;
537  deltaPrime.clear();
538  }
539  known.clear();
540  }
541 
542  if (oct->GetLocation() == eUnknown)
543  { // if the node does not know its location
544  vector<OctantSharedPtr> known;
545  map<OctantFace, vector<OctantSharedPtr> > nList =
546  oct->GetNeigbours();
547  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
548 
549  for (it = nList.begin(); it != nList.end(); it++)
550  {
551  for (int j = 0; j < it->second.size(); j++)
552  {
553  if (it->second[j]->GetLocation() != eUnknown)
554  {
555  known.push_back(it->second[j]);
556  }
557  }
558  }
559 
560  if (known.size() > 0)
561  {
562  vector<OctantSharedPtr> isNotOnBound;
563  for (int j = 0; j < known.size(); j++)
564  {
565  if (known[j]->GetLocation() != eOnBoundary)
566  {
567  isNotOnBound.push_back(known[j]);
568  }
569  }
570 
571  if (isNotOnBound.size() > 0)
572  {
573  oct->SetLocation(isNotOnBound[0]->GetLocation());
574  }
575  else
576  {
577  NekDouble dist = numeric_limits<double>::max();
578 
579  OctantSharedPtr closest;
580 
581  for (int j = 0; j < known.size(); j++)
582  {
583  if (oct->Distance(known[j]) < dist)
584  {
585  closest = known[j];
586  dist = oct->Distance(known[j]);
587  }
588  }
589 
590  CurvaturePointSharedPtr cp = closest->GetCPPoint();
591 
592  Array<OneD, NekDouble> octloc, cploc, vec(3), uv, N;
593  int surf;
594  cp->GetCAD(surf, uv);
595  N = m_cad->GetSurf(surf)->N(uv);
596 
597  octloc = oct->GetLoc();
598  cploc = cp->GetLoc();
599 
600  vec[0] = octloc[0] - cploc[0];
601  vec[1] = octloc[1] - cploc[1];
602  vec[2] = octloc[2] - cploc[2];
603 
604  NekDouble dot =
605  vec[0] * N[0] + vec[1] * N[1] + vec[2] * N[2];
606 
607  if (dot <= 0.0)
608  {
609  oct->SetLocation(eOutside);
610  ct += 1;
611  }
612  else
613  {
614  oct->SetLocation(eInside);
615  ct += 1;
616  }
617  }
618  }
619  known.clear();
620  }
621  }
622 
623  } while (ct > 0);
624 
625  for (int i = 0; i < m_octants.size(); i++)
626  {
627  ASSERTL0(m_octants[i]->IsDeltaKnown(),
628  "does not know delta after propergation");
629  }
630 }
631 
632 void Octree::SmoothAllOctants()
633 {
634  // until no more changes occur smooth the mesh specification between all
635  // m_octants not particualrly strictly
636  int ct = 0;
637 
638  do
639  {
640  ct = 0;
641  for (int i = 0; i < m_octants.size(); i++)
642  {
643  OctantSharedPtr oct = m_octants[i];
644 
645  vector<OctantSharedPtr> check;
646  map<OctantFace, vector<OctantSharedPtr> > nList =
647  oct->GetNeigbours();
648  map<OctantFace, vector<OctantSharedPtr> >::iterator it;
649  for (it = nList.begin(); it != nList.end(); it++)
650  {
651  for (int j = 0; j < it->second.size(); j++)
652  {
653  if (it->second[j]->GetDelta() < oct->GetDelta() &&
654  ddx(oct, it->second[j]) > 0.2)
655  {
656  check.push_back(it->second[j]);
657  }
658  }
659  }
660 
661  if (check.size() > 0)
662  {
663  NekDouble deltaSM = numeric_limits<double>::max();
664  for (int j = 0; j < check.size(); j++)
665  {
666  NekDouble r = oct->Distance(check[j]);
667 
668  if (0.199 * r + check[j]->GetDelta() < deltaSM)
669  {
670  deltaSM = 0.199 * r + check[j]->GetDelta();
671  }
672  }
673  oct->SetDelta(deltaSM);
674  ct += 1;
675  }
676  }
677 
678  } while (ct > 0);
679 }
680 
681 int Octree::CountElemt()
682 {
683  // by considering the volume of a tet evaluate the number of elements in the
684  // mesh
685 
686  NekDouble total = 0.0;
687 
688  Array<OneD, NekDouble> boundingBox = m_cad->GetBoundingBox();
689 
690  for (int i = 0; i < m_octants.size(); i++)
691  {
692  OctantSharedPtr oct = m_octants[i];
693  if (oct->GetLocation() == eInside)
694  {
695  total += 8.0 * oct->DX() * oct->DX() * oct->DX() /
696  (oct->GetDelta() * oct->GetDelta() * oct->GetDelta() /
697  6.0 / sqrt(2));
698  }
699  else if (oct->GetLocation() == eOnBoundary)
700  {
701  NekDouble vol = 1.0;
702  if (oct->FX(eBack) < boundingBox[1] &&
703  oct->FX(eForward) > boundingBox[0])
704  {
705  // then there is some over lap in x
706  NekDouble min, max;
707  if (oct->FX(eBack) > boundingBox[0])
708  {
709  min = oct->FX(eBack);
710  }
711  else
712  {
713  min = boundingBox[0];
714  }
715  if (boundingBox[1] < oct->FX(eForward))
716  {
717  max = boundingBox[1];
718  }
719  else
720  {
721  max = oct->FX(eForward);
722  }
723  vol *= (max - min);
724  }
725  else
726  {
727  vol *= 0.0;
728  }
729 
730  if (oct->FX(eDown) < boundingBox[3] &&
731  oct->FX(eUp) > boundingBox[2])
732  {
733  // then there is some over lap in x
734  NekDouble min, max;
735  if (oct->FX(eDown) > boundingBox[2])
736  {
737  min = oct->FX(eDown);
738  }
739  else
740  {
741  min = boundingBox[2];
742  }
743  if (boundingBox[3] < oct->FX(eUp))
744  {
745  max = boundingBox[3];
746  }
747  else
748  {
749  max = oct->FX(eUp);
750  }
751  vol *= (max - min);
752  }
753  else
754  {
755  vol *= 0.0;
756  }
757 
758  if (oct->FX(eRight) < boundingBox[5] &&
759  oct->FX(eLeft) > boundingBox[4])
760  {
761  // then there is some over lap in x
762  NekDouble min, max;
763  if (oct->FX(eRight) > boundingBox[4])
764  {
765  min = oct->FX(eRight);
766  }
767  else
768  {
769  min = boundingBox[4];
770  }
771  if (boundingBox[5] < oct->FX(eLeft))
772  {
773  max = boundingBox[5];
774  }
775  else
776  {
777  max = oct->FX(eLeft);
778  }
779  vol *= (max - min);
780  }
781  else
782  {
783  vol *= 0.0;
784  }
785  total += vol / 2.0 / (oct->GetDelta() * oct->GetDelta() *
786  oct->GetDelta() / 6.0 / sqrt(2));
787  }
788  }
789 
790  return int(total);
791 }
792 
794 {
796  NekDouble R, delta;
799  NekDouble r,
800  NekDouble d)
801  : x1(p1), x2(p2), R(r), delta(d)
802  {
803  }
804 
806  {
807  Array<OneD, NekDouble> Le(3), Re(3), s(3);
808  for (int i = 0; i < 3; i++)
809  {
810  Le[i] = p[i] - x1[i];
811  Re[i] = p[i] - x2[i];
812  s[i] = x2[i] - x1[i];
813  }
814  Array<OneD, NekDouble> dev(3);
815  dev[0] = Le[1] * Re[2] - Re[1] * Le[2];
816  dev[1] = Le[0] * Re[2] - Re[0] * Le[2];
817  dev[2] = Le[0] * Re[1] - Re[0] * Le[1];
818 
819  NekDouble dist =
820  sqrt(dev[0] * dev[0] + dev[1] * dev[1] + dev[2] * dev[2]) /
821  sqrt(s[0] * s[0] + s[1] * s[1] + s[2] * s[2]);
822 
823  NekDouble t = -1.0 * ((x1[0] - p[0]) * s[0] + (x1[1] - p[1]) * s[1] +
824  (x1[1] - p[1]) * s[1]) /
825  Length() / Length();
826 
827  if (dist < R && !(t > 1) && !(t < 0))
828  {
829  return true;
830  }
831  else
832  {
833  return false;
834  }
835  }
836 
838  {
839  return sqrt((x1[0] - x2[0]) * (x1[0] - x2[0]) +
840  (x1[1] - x2[1]) * (x1[1] - x2[1]) +
841  (x1[2] - x2[2]) * (x1[2] - x2[2]));
842  }
843 };
844 
845 void Octree::CompileCuravturePointList()
846 {
847  for (int i = 1; i <= m_cad->GetNumSurf(); i++)
848  {
849  CADSurfSharedPtr surf = m_cad->GetSurf(i);
850  Array<OneD, NekDouble> bounds = surf->GetBounds();
851 
852  // to figure out the amount of curvature sampling to be conducted on
853  // each parameter plane the surface is first sampled with a 40x40 grid
854  // the real space lengths of this grid are analysed to find the largest
855  // strecthing in the u and v directions
856  // this stretching this then cosnidered with the mindelta user input
857  // to find a number of sampling points in each direction which
858  // enures that in the final octree each surface octant will have at
859  // least
860  // one sample point within its volume.
861  // the 40x40 grid is used to ensure each surface has a minimum of 40x40
862  // samples.
863  NekDouble du = (bounds[1] - bounds[0]) / (40 - 1);
864  NekDouble dv = (bounds[3] - bounds[2]) / (40 - 1);
865 
866  NekDouble DeltaU = 0.0;
867  NekDouble DeltaV = 0.0;
868 
869  Array<TwoD, Array<OneD, NekDouble> > samplepoints(40, 40);
870 
871  for (int j = 0; j < 40; j++)
872  {
873  for (int k = 0; k < 40; k++)
874  {
876  uv[0] = k * du + bounds[0];
877  uv[1] = j * dv + bounds[2];
878  if (j == 40 - 1)
879  uv[1] = bounds[3];
880  if (k == 40 - 1)
881  uv[0] = bounds[1];
882  samplepoints[k][j] = surf->P(uv);
883  }
884  }
885 
886  for (int j = 0; j < 40 - 1; j++)
887  {
888  for (int k = 0; k < 40 - 1; k++)
889  {
890  NekDouble deltau = sqrt(
891  (samplepoints[k][j][0] - samplepoints[k + 1][j][0]) *
892  (samplepoints[k][j][0] - samplepoints[k + 1][j][0]) +
893  (samplepoints[k][j][1] - samplepoints[k + 1][j][1]) *
894  (samplepoints[k][j][1] - samplepoints[k + 1][j][1]) +
895  (samplepoints[k][j][2] - samplepoints[k + 1][j][2]) *
896  (samplepoints[k][j][2] - samplepoints[k + 1][j][2]));
897  NekDouble deltav = sqrt(
898  (samplepoints[k][j][0] - samplepoints[k][j + 1][0]) *
899  (samplepoints[k][j][0] - samplepoints[k][j + 1][0]) +
900  (samplepoints[k][j][1] - samplepoints[k][j + 1][1]) *
901  (samplepoints[k][j][1] - samplepoints[k][j + 1][1]) +
902  (samplepoints[k][j][2] - samplepoints[k][j + 1][2]) *
903  (samplepoints[k][j][2] - samplepoints[k][j + 1][2]));
904 
905  if (deltau > DeltaU)
906  DeltaU = deltau;
907  if (deltav > DeltaV)
908  DeltaV = deltav;
909  }
910  }
911 
912  // these are the acutal number of sample points in each parametric
913  // direction
914  int nu = ceil(DeltaU / m_minDelta) * 40;
915  int nv = ceil(DeltaV / m_minDelta) * 40;
916 
917  for (int j = 0; j < nu; j++)
918  {
919  for (int k = 0; k < nv; k++)
920  {
922  uv[0] = (bounds[1] - bounds[0]) / (nu - 1) * j + bounds[0];
923  uv[1] = (bounds[3] - bounds[2]) / (nv - 1) * k + bounds[2];
924 
925  // this prevents round off error at the end of the surface
926  // may not be neseercary but works
927  if (j == nu - 1)
928  uv[0] = bounds[1];
929  if (k == nv - 1)
930  uv[1] = bounds[3];
931 
932  NekDouble C = surf->Curvature(uv);
933 
934  // create new point based on smallest R, flat surfaces have k=0
935  // but still need a point for element estimation
936  if (C != 0.0)
937  {
938  bool minlimited = false;
939  NekDouble ideal;
940 
941  NekDouble del =
942  2.0 * (1.0 / C) * sqrt(m_eps * (2.0 - m_eps));
943 
944  if (del > m_maxDelta)
945  {
946  del = m_maxDelta;
947  }
948  if (del < m_minDelta)
949  {
950  ideal = del;
951  del = m_minDelta;
952  minlimited = true;
953  }
954 
955  if (minlimited)
956  {
957  CurvaturePointSharedPtr newCPoint =
959  surf->GetId(), uv, surf->P(uv), del, ideal);
960 
961  m_cpList.push_back(newCPoint);
962  }
963  else
964  {
965  CurvaturePointSharedPtr newCPoint =
967  surf->GetId(), uv, surf->P(uv), del);
968 
969  m_cpList.push_back(newCPoint);
970  }
971  }
972  else
973  {
974  CurvaturePointSharedPtr newCPoint =
976  surf->GetId(), uv, surf->P(uv));
977 
978  m_cpList.push_back(newCPoint);
979  }
980  }
981  }
982  }
983 
984  if (m_udsfile == "N")
985  {
986  return;
987  }
988 
989  // now deal with the user defined spacing
990  vector<linesource> lsources;
991  fstream fle;
992  fle.open(m_udsfile.c_str());
993 
994  string fileline;
995 
996  while (!fle.eof())
997  {
998  getline(fle, fileline);
999  stringstream s(fileline);
1000  string word;
1001  s >> word;
1002  if (word == "#")
1003  {
1004  continue;
1005  }
1006 
1007  Array<OneD, NekDouble> x1(3), x2(3);
1008  NekDouble r, d;
1009  x1[0] = boost::lexical_cast<double>(word);
1010  s >> x1[1] >> x1[2] >> x2[0] >> x2[1] >> x2[2] >> r >> d;
1011 
1012  lsources.push_back(linesource(x1, x2, r, d));
1013  }
1014  fle.close();
1015 
1016  for (int j = 0; j < lsources.size(); j++)
1017  {
1018  cout << lsources[j].x1[0] << " " << lsources[j].x1[1] << " "
1019  << lsources[j].x1[2] << endl;
1020  cout << lsources[j].x2[0] << " " << lsources[j].x2[1] << " "
1021  << lsources[j].x2[2] << endl;
1022  cout << lsources[j].Length() << endl;
1023  }
1024 
1025  int ct = 0;
1026  for (int i = 0; i < m_cpList.size(); i++)
1027  {
1028  for (int j = 0; j < lsources.size(); j++)
1029  {
1030  if (lsources[j].withinRange(m_cpList[i]->GetLoc()))
1031  {
1032  ct++;
1033  m_cpList[i]->SetDelta(lsources[j].delta);
1034  }
1035  }
1036  }
1037  cout << ct << endl;
1038 
1039  ///@TODO need to add curvature points with the false tag to make octree
1040  ///modification work
1041  // off surfaces
1042 
1043  /*for(int i = 0; i < lsources.size(); i++)
1044  {
1045  int nc; //number of point to add cicularly
1046  int nl; //number of point to add length
1047 
1048  nc = ceil(2.0*3.142*lsources[i].R / lsources[i].delta)*2;
1049  nl = ceil(lsources[i].Length() / lsources[i].delta)*2;
1050 
1051  NekDouble dr = lsources[i].Length() / nl;
1052  NekDouble dtheta = 2.0*3.142 / nc;
1053 
1054  for(int j = 0; j < nl; j++)
1055  {
1056  NekDouble len =
1057  for(int k = 0; k < nc; k++)
1058  {
1059 
1060  }
1061  }
1062 
1063  }*/
1064 }
1065 
1067 {
1068  return fabs(i->GetDelta() - j->GetDelta()) / i->Distance(j);
1069 }
1070 }
1071 }
#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
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
STL namespace.
Array< OneD, NekDouble > x2
Definition: Octree.cpp:795
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
Represents a point in the domain.
Definition: Node.h:60
boost::shared_ptr< CurvaturePoint > CurvaturePointSharedPtr
double NekDouble
boost::shared_ptr< Octant > OctantSharedPtr
Definition: Octant.h:74
linesource(Array< OneD, NekDouble > p1, Array< OneD, NekDouble > p2, NekDouble r, NekDouble d)
Definition: Octree.cpp:797
boost::shared_ptr< CADSurf > CADSurfSharedPtr
Definition: CADSurf.h:217
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:137
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:52
bool withinRange(Array< OneD, NekDouble > p)
Definition: Octree.cpp:805