Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Octant.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Octant.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: octant object methods.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
37 
38 using namespace std;
39 namespace Nektar
40 {
41 namespace NekMeshUtils
42 {
43 
45 {
46  switch (f)
47  {
48  case eUp:
49  return eDown;
50  case eDown:
51  return eUp;
52  case eForward:
53  return eBack;
54  case eBack:
55  return eForward;
56  case eLeft:
57  return eRight;
58  case eRight:
59  return eLeft;
60  }
61 
62  return eUp;
63 }
64 
65 Octant::Octant(int i, OctantSharedPtr p, Array<OneD, OctantFace> dir)
66  : m_id(i), m_parent(p)
67 {
68  // initialise variables to defualt states
69  m_leaf = true;
70  m_needToDivide = false;
71  m_numValidPoints = 0;
73  m_delta = pair<bool, NekDouble>(false, 0.0);
74  NekDouble maxDif = 0;
75  NekDouble minDif = numeric_limits<double>::max();
77 
78  // build empty neigbour map
79  m_neigbours[eUp] = vector<OctantSharedPtr>();
80  m_neigbours[eDown] = vector<OctantSharedPtr>();
81  m_neigbours[eForward] = vector<OctantSharedPtr>();
82  m_neigbours[eBack] = vector<OctantSharedPtr>();
83  m_neigbours[eLeft] = vector<OctantSharedPtr>();
84  m_neigbours[eRight] = vector<OctantSharedPtr>();
85 
86  // pull information from parent
87  Array<OneD, NekDouble> parentloc = m_parent->GetLoc();
88  m_hd = m_parent->DX() / 2.0;
90  if (dir[0] == eForward)
91  {
92  m_loc[0] = parentloc[0] + m_hd;
93  }
94  else
95  {
96  m_loc[0] = parentloc[0] - m_hd;
97  }
98  if (dir[1] == eUp)
99  {
100  m_loc[1] = parentloc[1] + m_hd;
101  }
102  else
103  {
104  m_loc[1] = parentloc[1] - m_hd;
105  }
106  if (dir[2] == eLeft)
107  {
108  m_loc[2] = parentloc[2] + m_hd;
109  }
110  else
111  {
112  m_loc[2] = parentloc[2] - m_hd;
113  }
114 
115  vector<SPBaseSharedPtr> SourcePointList = m_parent->GetSPList();
116 
117  // setup complete
118 
119  // look over the curvature point list provided by the parent,
120  // firstly look to see if it is in the new octant and if so
121  // add it to the conserdation of the delta specification
122  for (int i = 0; i < SourcePointList.size(); i++)
123  {
124  Array<OneD, NekDouble> cploc = SourcePointList[i]->GetLoc();
125  if (!(cploc[0] > m_loc[0] + m_hd || cploc[0] < m_loc[0] - m_hd ||
126  cploc[1] > m_loc[1] + m_hd || cploc[1] < m_loc[1] - m_hd ||
127  cploc[2] > m_loc[2] + m_hd || cploc[2] < m_loc[2] - m_hd))
128  {
129  m_localSPList.push_back(SourcePointList[i]);
130 
131  if (SourcePointList[i]->HasDelta())
132  {
133  if (SourcePointList[i]->GetDelta() > maxDif)
134  {
135  maxDif = SourcePointList[i]->GetDelta();
136  }
137 
138  if (SourcePointList[i]->GetDelta() < minDif)
139  {
140  minDif = SourcePointList[i]->GetDelta();
141  }
143  }
144  if (SourcePointList[i]->Isboundary())
145  {
147  }
148  }
149  }
150 
151  // if it has valid points delta can be calculated
152  if (NumValidCurvePoint() > 0)
153  {
154  // geometrically octant should subdivide
155  if (maxDif / minDif > 1.5)
156  {
157  m_needToDivide = true;
158  }
159 
160  SetDelta(minDif);
161 
162  //encourage subdivision to keep spec smooth
163  if (GetDelta() < 5.0 * DX())
164  {
165  m_needToDivide = true;
166  }
167  }
168 
169  if (GetNumBoundary() > 0)
170  {
172  }
173 }
174 
175 // constructor for the master octant
177  NekDouble x,
178  NekDouble y,
179  NekDouble z,
180  NekDouble dx,
181  const vector<SPBaseSharedPtr> &splist)
182  : m_id(i), m_hd(dx)
183 {
184  m_neigbours[eUp] = vector<OctantSharedPtr>();
185  m_neigbours[eDown] = vector<OctantSharedPtr>();
186  m_neigbours[eForward] = vector<OctantSharedPtr>();
187  m_neigbours[eBack] = vector<OctantSharedPtr>();
188  m_neigbours[eLeft] = vector<OctantSharedPtr>();
189  m_neigbours[eRight] = vector<OctantSharedPtr>();
190 
191  // initialise variables to defualt states
192  m_leaf = true;
193  m_needToDivide = true;
194  m_numValidPoints = 0;
195  m_delta = pair<bool, NekDouble>(false, 0.0);
196 
198  m_loc[0] = x;
199  m_loc[1] = y;
200  m_loc[2] = z;
201 
202  m_localSPList = splist;
203 
204  for (int i = 0; i < m_localSPList.size(); i++)
205  {
206  if (m_localSPList[i]->HasDelta())
207  {
209  }
210  }
211 
213 }
214 
216 {
217  ASSERTL0(m_leaf, "octant must be a leaf for subdivision");
218 
219  m_leaf = false; // set as not leaf and make children
220 
221  // need to loop over all neigbours and remove this octant from their lists
222  for (int i = 0; i < 6; i++)
223  {
224  OctantFace f = static_cast<OctantFace>(i);
225  vector<OctantSharedPtr> os = m_neigbours[f];
226  for (int j = 0; j < os.size(); j++)
227  {
228  os[j]->RemoveNeigbour(GetId(), GetReverseFace(f));
229  }
230  }
231 
232  Array<OneD, OctantSharedPtr> children(8);
233 
234  for (int i = 0; i < 8; i++)
235  {
236  // set up x,y,z ordering of the 8 octants
238  if (i < 4)
239  {
240  dir[0] = eForward;
241  if (i < 2)
242  {
243  dir[1] = eUp;
244  }
245  else
246  {
247  dir[1] = eDown;
248  }
249  if (i == 0 || i == 2)
250  {
251  dir[2] = eLeft;
252  }
253  else
254  {
255  dir[2] = eRight;
256  }
257  }
258  else
259  {
260  dir[0] = eBack;
261  if (i < 6)
262  {
263  dir[1] = eUp;
264  }
265  else
266  {
267  dir[1] = eDown;
268  }
269  if (i == 4 || i == 6)
270  {
271  dir[2] = eLeft;
272  }
273  else
274  {
275  dir[2] = eRight;
276  }
277  }
278 
279  children[i] = boost::shared_ptr<Octant>(new Octant(numoct++, p, dir));
280  }
281 
282  SetChildren(children);
283 
284  // this set of neibours are based on the children of the octant, only covers
285  // three sides
286  children[0]->SetNeigbour(children[1], eRight);
287  children[0]->SetNeigbour(children[4], eBack);
288  children[0]->SetNeigbour(children[2], eDown);
289 
290  children[1]->SetNeigbour(children[0], eLeft);
291  children[1]->SetNeigbour(children[5], eBack);
292  children[1]->SetNeigbour(children[3], eDown);
293 
294  children[2]->SetNeigbour(children[3], eRight);
295  children[2]->SetNeigbour(children[6], eBack);
296  children[2]->SetNeigbour(children[0], eUp);
297 
298  children[3]->SetNeigbour(children[2], eLeft);
299  children[3]->SetNeigbour(children[7], eBack);
300  children[3]->SetNeigbour(children[1], eUp);
301 
302  children[4]->SetNeigbour(children[5], eRight);
303  children[4]->SetNeigbour(children[0], eForward);
304  children[4]->SetNeigbour(children[6], eDown);
305 
306  children[5]->SetNeigbour(children[4], eLeft);
307  children[5]->SetNeigbour(children[1], eForward);
308  children[5]->SetNeigbour(children[7], eDown);
309 
310  children[6]->SetNeigbour(children[7], eRight);
311  children[6]->SetNeigbour(children[2], eForward);
312  children[6]->SetNeigbour(children[4], eUp);
313 
314  children[7]->SetNeigbour(children[6], eLeft);
315  children[7]->SetNeigbour(children[3], eForward);
316  children[7]->SetNeigbour(children[5], eUp);
317 
318  // need to obtain the remaning information from the parents neigbours
319  // (m_neigbours)
320  // consider top face
321  if (m_neigbours[eUp].size() == 1)
322  {
323  children[0]->SetNeigbour(m_neigbours[eUp][0], eUp);
324  children[1]->SetNeigbour(m_neigbours[eUp][0], eUp);
325  children[4]->SetNeigbour(m_neigbours[eUp][0], eUp);
326  children[5]->SetNeigbour(m_neigbours[eUp][0], eUp);
327  m_neigbours[eUp][0]->SetNeigbour(children[0], eDown);
328  m_neigbours[eUp][0]->SetNeigbour(children[1], eDown);
329  m_neigbours[eUp][0]->SetNeigbour(children[4], eDown);
330  m_neigbours[eUp][0]->SetNeigbour(children[5], eDown);
331  }
332  else if (m_neigbours[eUp].size() == 4)
333  {
334  children[0]->SetNeigbour(m_neigbours[eUp][0], eUp); // 2
335  children[1]->SetNeigbour(m_neigbours[eUp][1], eUp); // 3
336  children[4]->SetNeigbour(m_neigbours[eUp][2], eUp); // 6
337  children[5]->SetNeigbour(m_neigbours[eUp][3], eUp); // 7
338  m_neigbours[eUp][0]->SetNeigbour(children[0], eDown);
339  m_neigbours[eUp][1]->SetNeigbour(children[1], eDown);
340  m_neigbours[eUp][2]->SetNeigbour(children[4], eDown);
341  m_neigbours[eUp][3]->SetNeigbour(children[5], eDown);
342  }
343  else if (m_neigbours[eUp].size() != 0)
344  {
345  cout << "!!!!!"
346  << "NOT GOOD"
347  << "!!!!! " << m_neigbours[eUp].size() << endl;
348  }
349 
350  if (m_neigbours[eDown].size() == 1)
351  {
352  children[2]->SetNeigbour(m_neigbours[eDown][0], eDown);
353  children[3]->SetNeigbour(m_neigbours[eDown][0], eDown);
354  children[6]->SetNeigbour(m_neigbours[eDown][0], eDown);
355  children[7]->SetNeigbour(m_neigbours[eDown][0], eDown);
356  m_neigbours[eDown][0]->SetNeigbour(children[2], eUp);
357  m_neigbours[eDown][0]->SetNeigbour(children[3], eUp);
358  m_neigbours[eDown][0]->SetNeigbour(children[6], eUp);
359  m_neigbours[eDown][0]->SetNeigbour(children[7], eUp);
360  }
361  else if (m_neigbours[eDown].size() == 4)
362  {
363  children[2]->SetNeigbour(m_neigbours[eDown][0], eDown); // 0
364  children[3]->SetNeigbour(m_neigbours[eDown][1], eDown); // 1
365  children[6]->SetNeigbour(m_neigbours[eDown][2], eDown); // 4
366  children[7]->SetNeigbour(m_neigbours[eDown][3], eDown); // 5
367  m_neigbours[eDown][0]->SetNeigbour(children[2], eUp);
368  m_neigbours[eDown][1]->SetNeigbour(children[3], eUp);
369  m_neigbours[eDown][2]->SetNeigbour(children[6], eUp);
370  m_neigbours[eDown][3]->SetNeigbour(children[7], eUp);
371  }
372  else if (m_neigbours[eDown].size() != 0)
373  {
374  cout << "!!!!!"
375  << "NOT GOOD"
376  << "!!!!! " << m_neigbours[eDown].size() << endl;
377  }
378 
379  if (m_neigbours[eForward].size() == 1)
380  {
381  children[0]->SetNeigbour(m_neigbours[eForward][0], eForward);
382  children[1]->SetNeigbour(m_neigbours[eForward][0], eForward);
383  children[2]->SetNeigbour(m_neigbours[eForward][0], eForward);
384  children[3]->SetNeigbour(m_neigbours[eForward][0], eForward);
385  m_neigbours[eForward][0]->SetNeigbour(children[0], eBack);
386  m_neigbours[eForward][0]->SetNeigbour(children[1], eBack);
387  m_neigbours[eForward][0]->SetNeigbour(children[2], eBack);
388  m_neigbours[eForward][0]->SetNeigbour(children[3], eBack);
389  }
390  else if (m_neigbours[eForward].size() == 4)
391  {
392  children[0]->SetNeigbour(m_neigbours[eForward][0], eForward); // 4
393  children[1]->SetNeigbour(m_neigbours[eForward][1], eForward); // 5
394  children[2]->SetNeigbour(m_neigbours[eForward][2], eForward); // 6
395  children[3]->SetNeigbour(m_neigbours[eForward][3], eForward); // 7
396  m_neigbours[eForward][0]->SetNeigbour(children[0], eBack);
397  m_neigbours[eForward][1]->SetNeigbour(children[1], eBack);
398  m_neigbours[eForward][2]->SetNeigbour(children[2], eBack);
399  m_neigbours[eForward][3]->SetNeigbour(children[3], eBack);
400  }
401  else if (m_neigbours[eForward].size() != 0)
402  {
403  cout << "!!!!!"
404  << "NOT GOOD"
405  << "!!!!! " << m_neigbours[eForward].size() << endl;
406  }
407 
408  if (m_neigbours[eBack].size() == 1)
409  {
410  children[4]->SetNeigbour(m_neigbours[eBack][0], eBack);
411  children[5]->SetNeigbour(m_neigbours[eBack][0], eBack);
412  children[6]->SetNeigbour(m_neigbours[eBack][0], eBack);
413  children[7]->SetNeigbour(m_neigbours[eBack][0], eBack);
414  m_neigbours[eBack][0]->SetNeigbour(children[4], eForward);
415  m_neigbours[eBack][0]->SetNeigbour(children[5], eForward);
416  m_neigbours[eBack][0]->SetNeigbour(children[6], eForward);
417  m_neigbours[eBack][0]->SetNeigbour(children[7], eForward);
418  }
419  else if (m_neigbours[eBack].size() == 4)
420  {
421  children[4]->SetNeigbour(m_neigbours[eBack][0], eBack); // 0
422  children[5]->SetNeigbour(m_neigbours[eBack][1], eBack); // 1
423  children[6]->SetNeigbour(m_neigbours[eBack][2], eBack); // 2
424  children[7]->SetNeigbour(m_neigbours[eBack][3], eBack); // 3
425  m_neigbours[eBack][0]->SetNeigbour(children[4], eForward);
426  m_neigbours[eBack][1]->SetNeigbour(children[5], eForward);
427  m_neigbours[eBack][2]->SetNeigbour(children[6], eForward);
428  m_neigbours[eBack][3]->SetNeigbour(children[7], eForward);
429  }
430  else if (m_neigbours[eBack].size() != 0)
431  {
432  cout << "!!!!!"
433  << "NOT GOOD"
434  << "!!!!! " << m_neigbours[eBack].size() << endl;
435  }
436 
437  if (m_neigbours[eLeft].size() == 1)
438  {
439  children[0]->SetNeigbour(m_neigbours[eLeft][0], eLeft);
440  children[2]->SetNeigbour(m_neigbours[eLeft][0], eLeft);
441  children[4]->SetNeigbour(m_neigbours[eLeft][0], eLeft);
442  children[6]->SetNeigbour(m_neigbours[eLeft][0], eLeft);
443  m_neigbours[eLeft][0]->SetNeigbour(children[0], eRight);
444  m_neigbours[eLeft][0]->SetNeigbour(children[2], eRight);
445  m_neigbours[eLeft][0]->SetNeigbour(children[4], eRight);
446  m_neigbours[eLeft][0]->SetNeigbour(children[6], eRight);
447  }
448  else if (m_neigbours[eLeft].size() == 4)
449  {
450  children[0]->SetNeigbour(m_neigbours[eLeft][0], eLeft); // 1
451  children[2]->SetNeigbour(m_neigbours[eLeft][1], eLeft); // 3
452  children[4]->SetNeigbour(m_neigbours[eLeft][2], eLeft); // 5
453  children[6]->SetNeigbour(m_neigbours[eLeft][3], eLeft); // 7
454  m_neigbours[eLeft][0]->SetNeigbour(children[0], eRight);
455  m_neigbours[eLeft][1]->SetNeigbour(children[2], eRight);
456  m_neigbours[eLeft][2]->SetNeigbour(children[4], eRight);
457  m_neigbours[eLeft][3]->SetNeigbour(children[6], eRight);
458  }
459  else if (m_neigbours[eLeft].size() != 0)
460  {
461  cout << "!!!!!"
462  << "NOT GOOD"
463  << "!!!!! " << m_neigbours[eLeft].size() << endl;
464  cout << m_neigbours[eLeft].size() << endl;
465  }
466 
467  if (m_neigbours[eRight].size() == 1)
468  {
469  children[1]->SetNeigbour(m_neigbours[eRight][0], eRight);
470  children[3]->SetNeigbour(m_neigbours[eRight][0], eRight);
471  children[5]->SetNeigbour(m_neigbours[eRight][0], eRight);
472  children[7]->SetNeigbour(m_neigbours[eRight][0], eRight);
473  m_neigbours[eRight][0]->SetNeigbour(children[1], eLeft);
474  m_neigbours[eRight][0]->SetNeigbour(children[3], eLeft);
475  m_neigbours[eRight][0]->SetNeigbour(children[5], eLeft);
476  m_neigbours[eRight][0]->SetNeigbour(children[7], eLeft);
477  }
478  else if (m_neigbours[eRight].size() == 4)
479  {
480  children[1]->SetNeigbour(m_neigbours[eRight][0], eRight); // 0
481  children[3]->SetNeigbour(m_neigbours[eRight][1], eRight); // 2
482  children[5]->SetNeigbour(m_neigbours[eRight][2], eRight); // 4
483  children[7]->SetNeigbour(m_neigbours[eRight][3], eRight); // 6
484  m_neigbours[eRight][0]->SetNeigbour(children[1], eLeft);
485  m_neigbours[eRight][1]->SetNeigbour(children[3], eLeft);
486  m_neigbours[eRight][2]->SetNeigbour(children[5], eLeft);
487  m_neigbours[eRight][3]->SetNeigbour(children[7], eLeft);
488  }
489  else if (m_neigbours[eRight].size() != 0)
490  {
491  cout << "!!!!!"
492  << "NOT GOOD"
493  << "!!!!! " << m_neigbours[eRight].size() << endl;
494  }
495 }
496 
498 {
499  vector<OctantSharedPtr> tmp = m_neigbours[f];
500  m_neigbours[f].clear();
501  for (int i = 0; i < tmp.size(); i++)
502  {
503  if (tmp[i]->GetId() != id)
504  {
505  m_neigbours[f].push_back(tmp[i]);
506  }
507  }
508 }
509 
510 bool operator==(OctantSharedPtr const &p1, OctantSharedPtr const &p2)
511 {
512  if (p1->GetId() == p2->GetId())
513  {
514  return true;
515  }
516 
517  return false;
518 }
519 }
520 }
NekDouble GetDelta()
Get value of delta.
Definition: Octant.h:201
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
std::vector< SPBaseSharedPtr > m_localSPList
curvature sampling point list
Definition: Octant.h:344
OctantLocation m_location
idenify if delta has ben set
Definition: Octant.h:353
void RemoveNeigbour(int id, OctantFace f)
Remove a neigbour from this octants list.
Definition: Octant.cpp:497
bool m_leaf
leaf identifer
Definition: Octant.h:334
STL namespace.
OctantSharedPtr m_parent
parent id
Definition: Octant.h:336
void SetDelta(NekDouble d)
Set the value for delta for this octant.
Definition: Octant.h:192
NekDouble m_hd
half dimension of the octant
Definition: Octant.h:342
NekDouble DX()
Get the octants half dimension.
Definition: Octant.h:142
Octant(int i, OctantSharedPtr p, Array< OneD, OctantFace > dir)
Defualt constructor.
Definition: Octant.cpp:65
bool operator==(ElmtConfig const &c1, ElmtConfig const &c2)
Compares two element config structs.
void Subdivide(OctantSharedPtr p, int &numoct)
Subdivide the octant.
Definition: Octant.cpp:215
int m_numValidPoints
number of valid cp points
Definition: Octant.h:346
int GetId()
Get the Id of the octant.
Definition: Octant.h:126
int NumValidCurvePoint()
Get the number of valid cp points in the octants list.
Definition: Octant.h:184
OctantFace GetReverseFace(OctantFace f)
Definition: Octant.cpp:44
double NekDouble
Array< OneD, NekDouble > m_loc
x,y,z location of the octant
Definition: Octant.h:340
boost::shared_ptr< Octant > OctantSharedPtr
Definition: Octant.h:74
OctantFace
enumeration of the 6 faces of a cube/octant
Definition: Octant.h:52
bool m_needToDivide
idenify if division is needed
Definition: Octant.h:351
void SetChildren(Array< OneD, OctantSharedPtr > c)
Set the children of this octant.
Definition: Octant.h:213
std::map< OctantFace, std::vector< OctantSharedPtr > > m_neigbours
list of neighbours
Definition: Octant.h:355
std::pair< bool, NekDouble > m_delta
mesh sizing parameter
Definition: Octant.h:349