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