Octree.cpp 35.8 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
////////////////////////////////////////////////////////////////////////////////
//
//  File: octree.h
//
//  For more information, please see: http://www.nektar.info/
//
//  The MIT License
//
//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
//  Department of Aeronautics, Imperial College London (UK), and Scientific
//  Computing and Imaging Institute, University of Utah (USA).
//
//  License for the specific language governing rights and limitations under
//  Permission is hereby granted, free of charge, to any person obtaining a
//  copy of this software and associated documentation files (the "Software"),
//  to deal in the Software without restriction, including without limitation
//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
//  and/or sell copies of the Software, and to permit persons to whom the
//  Software is furnished to do so, subject to the following conditions:
//
//  The above copyright notice and this permission notice shall be included
//  in all copies or substantial portions of the Software.
//
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
//  DEALINGS IN THE SOFTWARE.
//
//  Description: cad object methods.
//
////////////////////////////////////////////////////////////////////////////////

36
#include "Octree.h"
37
#include <NekMeshUtils/CADSystem/CADCurve.h>
mike turner's avatar
mike turner committed
38
#include <NekMeshUtils/CADSystem/CADSurf.h>
39
#include <NekMeshUtils/Module/Module.h>
40

41
#include <LibUtilities/BasicUtils/ParseUtils.h>
42 43
#include <LibUtilities/BasicUtils/Progressbar.hpp>

Michael Turner's avatar
Michael Turner committed
44 45
#include <boost/algorithm/string.hpp>

46
using namespace std;
47 48
namespace Nektar
{
49
namespace NekMeshUtils
50
{
51

52 53 54 55 56 57 58 59
void Octree::Process()
{
    Array<OneD, NekDouble> boundingBox = m_mesh->m_cad->GetBoundingBox();

    // build curvature samples
    CompileSourcePointList();

    if (m_mesh->m_verbose)
Dave Moxey's avatar
Dave Moxey committed
60
    {
61
        cout << "\tCurvature samples: " << m_SPList.size() << endl;
Dave Moxey's avatar
Dave Moxey committed
62
    }
63

Michael Turner's avatar
Michael Turner committed
64
    // make master octant based on the bounding box of the domain
65 66 67 68 69 70 71 72 73 74 75 76 77
    m_dim = max((boundingBox[1] - boundingBox[0]) / 2.0,
                (boundingBox[3] - boundingBox[2]) / 2.0);

    m_dim = max(m_dim, (boundingBox[5] - boundingBox[4]) / 2.0);

    m_centroid    = Array<OneD, NekDouble>(3);
    m_centroid[0] = (boundingBox[1] + boundingBox[0]) / 2.0;
    m_centroid[1] = (boundingBox[3] + boundingBox[2]) / 2.0;
    m_centroid[2] = (boundingBox[5] + boundingBox[4]) / 2.0;

    m_masteroct = MemoryManager<Octant>::AllocateSharedPtr(
        0, m_centroid[0], m_centroid[1], m_centroid[2], m_dim, m_SPList);

mike turner's avatar
mike turner committed
78
    // begin recersive subdivision
79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
    SubDivide();

    m_octants.clear();
    m_masteroct->CompileLeaves(m_octants);

    if (m_mesh->m_verbose)
    {
        cout << "\tOctants: " << m_octants.size() << endl;
    }

    SmoothSurfaceOctants();

    PropagateDomain();

    SmoothAllOctants();

    if (m_mesh->m_verbose)
    {
        int elem = CountElemt();
        cout << "\tPredicted mesh: " << elem << endl;
    }
100 101
}

102 103
NekDouble Octree::Query(Array<OneD, NekDouble> loc)
{
104
    // starting at master octant 0 move through succsesive m_octants which
Michael Turner's avatar
Michael Turner committed
105
    // contain the point loc until a leaf is found
mike turner's avatar
mike turner committed
106
    // first search through sourcepoints
Michael Turner's avatar
Michael Turner committed
107

Michael Turner's avatar
Michael Turner committed
108 109
    NekDouble tmp = numeric_limits<double>::max();

mike turner's avatar
mike turner committed
110
    for (int i = 0; i < m_lsources.size(); i++)
Michael Turner's avatar
Michael Turner committed
111
    {
mike turner's avatar
mike turner committed
112
        if (m_lsources[i].withinRange(loc))
Michael Turner's avatar
Michael Turner committed
113
        {
mike turner's avatar
mike turner committed
114
            tmp = min(m_lsources[i].delta, tmp);
Michael Turner's avatar
Michael Turner committed
115 116 117
        }
    }

118
    OctantSharedPtr n = m_masteroct;
119
    int quad;
Michael Turner's avatar
Michael Turner committed
120

121
    bool found = false;
Michael Turner's avatar
Michael Turner committed
122

123
    while (!found)
124
    {
125
        Array<OneD, NekDouble> octloc = n->GetLoc();
Michael Turner's avatar
Michael Turner committed
126

127 128
        if (!(loc[0] < octloc[0]) && // forward
            !(loc[1] < octloc[1]) && // forward
mike turner's avatar
mike turner committed
129
            !(loc[2] < octloc[2]))   // forward
Michael Turner's avatar
Michael Turner committed
130
        {
131
            quad = 0;
132
        }
133 134
        else if (!(loc[0] < octloc[0]) && // forward
                 !(loc[1] < octloc[1]) && // forward
mike turner's avatar
mike turner committed
135
                 !(loc[2] > octloc[2]))   // back
136
        {
137
            quad = 1;
138
        }
139 140
        else if (!(loc[0] < octloc[0]) && // forward
                 !(loc[1] > octloc[1]) && // back
mike turner's avatar
mike turner committed
141
                 !(loc[2] < octloc[2]))   // forward
142
        {
143
            quad = 2;
144
        }
145 146
        else if (!(loc[0] < octloc[0]) && // forward
                 !(loc[1] > octloc[1]) && // back
mike turner's avatar
mike turner committed
147
                 !(loc[2] > octloc[2]))   // back
148
        {
149
            quad = 3;
150
        }
151 152
        else if (!(loc[0] > octloc[0]) && // back
                 !(loc[1] < octloc[1]) && // forward
mike turner's avatar
mike turner committed
153
                 !(loc[2] < octloc[2]))   // forward
154
        {
155
            quad = 4;
156
        }
157 158
        else if (!(loc[0] > octloc[0]) && // back
                 !(loc[1] < octloc[1]) && // forward
mike turner's avatar
mike turner committed
159
                 !(loc[2] > octloc[2]))   // back
160
        {
161
            quad = 5;
162
        }
163 164
        else if (!(loc[0] > octloc[0]) && // back
                 !(loc[1] > octloc[1]) && // back
mike turner's avatar
mike turner committed
165
                 !(loc[2] < octloc[2]))   // forward
166
        {
167
            quad = 6;
168
        }
169 170
        else if (!(loc[0] > octloc[0]) && // back
                 !(loc[1] > octloc[1]) && // back
mike turner's avatar
mike turner committed
171
                 !(loc[2] > octloc[2]))   // back
172
        {
173
            quad = 7;
174
        }
Michael Turner's avatar
Michael Turner committed
175 176
        else
        {
177
            ASSERTL0(false, "Cannot locate quadrant");
Michael Turner's avatar
Michael Turner committed
178
        }
Michael Turner's avatar
Michael Turner committed
179

180
        n = n->GetChild(quad);
Michael Turner's avatar
Michael Turner committed
181

182
        if (n->IsLeaf())
183
        {
184
            found = true;
Michael Turner's avatar
Michael Turner committed
185 186
        }
    }
's avatar
fix  
committed
187

mike turner's avatar
mike turner committed
188
    return min(n->GetDelta(), tmp);
189
}
190

's avatar
fix  
committed
191 192 193 194
NekDouble Octree::GetMinDelta()
{
    NekDouble tmp = numeric_limits<double>::max();

mike turner's avatar
mike turner committed
195
    for (int i = 0; i < m_lsources.size(); i++)
's avatar
fix  
committed
196
    {
mike turner's avatar
mike turner committed
197
        tmp = min(m_lsources[i].delta, tmp);
's avatar
fix  
committed
198
    }
mike turner's avatar
mike turner committed
199
    return min(m_minDelta, tmp);
's avatar
fix  
committed
200 201
}

202
void Octree::WriteOctree(string nm)
203
{
204
    MeshSharedPtr oct = std::shared_ptr<Mesh>(new Mesh());
205 206 207 208
    oct->m_expDim     = 3;
    oct->m_spaceDim   = 3;
    oct->m_nummode    = 2;

209
    for (int i = 0; i < m_octants.size(); i++)
210
    {
211
        /*if(m_octants[i]->GetLocation() != eOnBoundary)
212 213 214
        {
            continue;
        }*/
215 216 217

        vector<NodeSharedPtr> ns(8);

218 219 220
        ns[0] = std::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eBack),
                                               m_octants[i]->FX(eDown),
                                               m_octants[i]->FX(eRight)));
221

222 223 224
        ns[1] = std::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eForward),
                                               m_octants[i]->FX(eDown),
                                               m_octants[i]->FX(eRight)));
225

226 227 228
        ns[2] = std::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eForward),
                                               m_octants[i]->FX(eUp),
                                               m_octants[i]->FX(eRight)));
229

230 231 232
        ns[3] = std::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eBack),
                                               m_octants[i]->FX(eUp),
                                               m_octants[i]->FX(eRight)));
233

234 235 236
        ns[4] = std::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eBack),
                                               m_octants[i]->FX(eDown),
                                               m_octants[i]->FX(eLeft)));
237

238 239 240
        ns[5] = std::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eForward),
                                               m_octants[i]->FX(eDown),
                                               m_octants[i]->FX(eLeft)));
241

242 243 244
        ns[6] = std::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eForward),
                                               m_octants[i]->FX(eUp),
                                               m_octants[i]->FX(eLeft)));
245

246 247 248
        ns[7] = std::shared_ptr<Node>(new Node(0, m_octants[i]->FX(eBack),
                                               m_octants[i]->FX(eUp),
                                               m_octants[i]->FX(eLeft)));
249 250 251

        vector<int> tags;
        tags.push_back(0);
252
        ElmtConfig conf(LibUtilities::eHexahedron, 1, false, false);
253
        ElementSharedPtr E = GetElementFactory().CreateInstance(
254
            LibUtilities::eHexahedron, conf, ns, tags);
255
        oct->m_element[3].push_back(E);
256
    }
Michael Turner's avatar
Michael Turner committed
257

mike turner's avatar
mike turner committed
258 259
    ModuleSharedPtr mod =
        GetModuleFactory().CreateInstance(ModuleKey(eOutputModule, "xml"), oct);
260 261 262 263 264 265 266
    mod->RegisterConfig("outfile", nm);
    mod->ProcessVertices();
    mod->ProcessEdges();
    mod->ProcessFaces();
    mod->ProcessElements();
    mod->ProcessComposites();
    mod->Process();
Michael Turner's avatar
Michael Turner committed
267
}
Michael Turner's avatar
Michael Turner committed
268

269
void Octree::SubDivide()
Michael Turner's avatar
Michael Turner committed
270
{
271
    bool repeat;
272 273
    int ct = 1;

274 275 276
    m_numoct = 1;
    m_masteroct->Subdivide(m_masteroct, m_numoct);

277
    if (m_mesh->m_verbose)
278 279 280
    {
        cout << "\tSubdivide iteration: ";
    }
Michael Turner's avatar
Michael Turner committed
281 282 283

    do
    {
284
        if (m_mesh->m_verbose)
285
        {
286 287 288
            cout << "\r                                                       ";
            cout << "\r";
            cout << "\tSubdivide iteration: " << ct;
289 290
            cout.flush();
        }
291
        ct++;
292
        repeat = false;
293
        m_octants.clear();
294
        // grab a list of the leaves curently in the octree
295
        m_masteroct->CompileLeaves(m_octants);
296

297
        VerifyNeigbours();
298

299 300 301 302 303 304 305 306
        // neeed to create a divide list, in the first list will be m_octants
        // which need to
        // subdivide based on curvature,
        // in the next list will be ocants which need to subdivide to make sure
        // level criteria is statisified for the previous list and so on
        // the list will keep building till no more lists are required.
        // the list will then be iterated through backwards to subdivide all the
        // sublists in turn.
307 308
        vector<vector<OctantSharedPtr> > dividelist;
        set<int> inlist;
309
        // build initial list
310
        {
311
            vector<OctantSharedPtr> sublist;
312
            for (int i = 0; i < m_octants.size(); i++)
313
            {
314 315
                if (m_octants[i]->NeedDivide() &&
                    m_octants[i]->DX() / 4.0 > m_minDelta)
316
                {
317 318
                    sublist.push_back(m_octants[i]);
                    inlist.insert(m_octants[i]->GetId());
319 320
                    repeat = true; // if there is a subdivision this whole
                                   // process needs to be repeated
321 322
                }
            }
323
            dividelist.push_back(sublist);
324
        }
325
        // then loop over building sublists until no more are required
326
        int ct2 = 0;
327
        while (true)
Michael Turner's avatar
Michael Turner committed
328
        {
329
            ct2++;
330 331 332
            vector<OctantSharedPtr> newsublist,
                previouslist = dividelist.back();
            for (int i = 0; i < previouslist.size(); i++)
Michael Turner's avatar
Michael Turner committed
333
            {
334 335
                map<OctantFace, vector<OctantSharedPtr> > nlist =
                    previouslist[i]->GetNeigbours();
336
                map<OctantFace, vector<OctantSharedPtr> >::iterator it;
337
                for (it = nlist.begin(); it != nlist.end(); it++)
Michael Turner's avatar
Michael Turner committed
338
                {
339
                    for (int j = 0; j < it->second.size(); j++)
Michael Turner's avatar
Michael Turner committed
340
                    {
341
                        if (previouslist[i]->DX() < it->second[j]->DX())
Michael Turner's avatar
Michael Turner committed
342
                        {
343 344 345
                            set<int>::iterator s =
                                inlist.find(it->second[j]->GetId());
                            if (s == inlist.end())
346 347
                            {
                                inlist.insert(it->second[j]->GetId());
348
                                newsublist.push_back(it->second[j]);
349
                            }
Michael Turner's avatar
Michael Turner committed
350 351 352 353
                        }
                    }
                }
            }
354
            if (newsublist.size() == 0)
355 356 357 358 359 360 361
            {
                break;
            }
            else
            {
                dividelist.push_back(newsublist);
            }
Michael Turner's avatar
Michael Turner committed
362
        }
363 364

        vector<vector<OctantSharedPtr> >::reverse_iterator rit;
365
        for (rit = dividelist.rbegin(); rit != dividelist.rend(); rit++)
366
        {
367
            vector<OctantSharedPtr> currentlist = *rit;
368
            for (int i = 0; i < currentlist.size(); i++)
369
            {
370
                currentlist[i]->Subdivide(currentlist[i], m_numoct);
371
            }
372
        }
373
    } while (repeat);
374

375
    if (m_mesh->m_verbose)
376 377 378
    {
        cout << endl;
    }
379
}
380

381
bool Octree::VerifyNeigbours()
382
{
Michael Turner's avatar
Michael Turner committed
383 384 385 386
    // check all octant links to their neighbours
    // at all times in the subdivision the set of neigbours must
    // conform to a set of criteria such as smoothness in size
    // this checks that
387
    bool error = false;
388
    for (int i = 0; i < m_octants.size(); i++)
389 390
    {
        bool valid = true;
391 392
        map<OctantFace, vector<OctantSharedPtr> > nlist =
            m_octants[i]->GetNeigbours();
393
        map<OctantFace, vector<OctantSharedPtr> >::iterator it;
394
        for (it = nlist.begin(); it != nlist.end(); it++)
395
        {
396
            if (it->second.size() == 0)
397
            {
mike turner's avatar
mike turner committed
398
                NekDouble expectedfx = 0.0;
399 400 401
                switch (it->first)
                {
                    case eUp:
402
                        expectedfx = m_centroid[1] + m_dim;
403 404
                        break;
                    case eDown:
405
                        expectedfx = m_centroid[1] - m_dim;
406 407
                        break;
                    case eLeft:
408
                        expectedfx = m_centroid[2] + m_dim;
409 410
                        break;
                    case eRight:
411
                        expectedfx = m_centroid[2] - m_dim;
412 413
                        break;
                    case eForward:
414
                        expectedfx = m_centroid[0] + m_dim;
415 416
                        break;
                    case eBack:
417
                        expectedfx = m_centroid[0] - m_dim;
418 419
                        break;
                }
420
                if (fabs(m_octants[i]->FX(it->first) - expectedfx) > 1E-6)
421 422 423
                {
                    valid = false;
                    cout << "wall neigbour error" << endl;
424 425
                    cout << expectedfx << " " << m_octants[i]->FX(it->first)
                         << " " << it->first << endl;
426 427
                }
            }
428
            else if (it->second.size() == 1)
429
            {
430 431
                if (!(m_octants[i]->DX() == it->second[0]->DX() ||
                      it->second[0]->DX() == 2.0 * m_octants[i]->DX()))
432 433 434
                {
                    valid = false;
                    cout << " 1 neigbour error" << endl;
435 436
                    cout << m_octants[i]->DX() << " " << it->second[0]->DX()
                         << endl;
437 438
                }
            }
439
            else if (it->second.size() == 4)
440
            {
441
                if (!(m_octants[i]->DX() / 2.0 == it->second[0]->DX()))
442 443 444
                {
                    valid = false;
                    cout << "4 neibour error" << endl;
445 446
                    cout << m_octants[i]->DX() << " " << it->second[0]->DX()
                         << endl;
447
                }
448
            }
449
        }
450
        if (!valid)
451 452 453 454
        {
            error = true;
            cout << "invalid neigbour config" << endl;
        }
455
    }
456
    return !error;
457 458
}

459
void Octree::SmoothSurfaceOctants()
460
{
461 462 463
    // for all the m_octants which are surface containing and know their delta
    // specification, look over all neighbours and ensure the specification
    // between them is smooth
464
    int ct = 0;
Michael Turner's avatar
Michael Turner committed
465

466 467
    do
    {
468 469
        ct = 0;
        for (int i = 0; i < m_octants.size(); i++)
470
        {
471
            OctantSharedPtr oct = m_octants[i];
472

473
            if (oct->IsDeltaKnown())
474
            {
475
                vector<OctantSharedPtr> check;
476
                map<OctantFace, vector<OctantSharedPtr> > nList =
477
                    oct->GetNeigbours();
478
                map<OctantFace, vector<OctantSharedPtr> >::iterator it;
Michael Turner's avatar
Michael Turner committed
479

480
                for (it = nList.begin(); it != nList.end(); it++)
481
                {
482
                    for (int j = 0; j < it->second.size(); j++)
483
                    {
484 485
                        if (it->second[j]->IsDeltaKnown() &&
                            it->second[j]->GetDelta() < oct->GetDelta() &&
Michael Turner's avatar
Michael Turner committed
486
                            ddx(oct, it->second[j]) > 0.2)
487 488 489
                        {
                            check.push_back(it->second[j]);
                        }
490
                    }
491
                }
Michael Turner's avatar
Michael Turner committed
492

493 494
                // for each neighbour listed in check_id, figure out the
                // smoothed
495
                // delta, and asign the miminum of these to nodes[i].GetDelta()
496
                if (check.size() > 0)
497 498
                {
                    NekDouble deltaSM = numeric_limits<double>::max();
499
                    for (int j = 0; j < check.size(); j++)
500
                    {
501
                        NekDouble r = oct->Distance(check[j]);
502

Michael Turner's avatar
Michael Turner committed
503
                        if (0.199 * r + check[j]->GetDelta() < deltaSM)
504
                        {
Michael Turner's avatar
Michael Turner committed
505
                            deltaSM = 0.199 * r + check[j]->GetDelta();
506 507
                        }
                    }
508
                    oct->SetDelta(deltaSM);
509
                    ct += 1;
510 511
                }
            }
512
        }
513
    } while (ct > 0);
514
}
Michael Turner's avatar
Michael Turner committed
515

516 517
void Octree::PropagateDomain()
{
518 519 520 521 522
    // until all m_octants know their delta specifcation and orientaion
    // look over all m_octants and if their neighours know either their
    // orientation
    // or specifcation calculate one for this octant
    int ct = 0;
Michael Turner's avatar
Michael Turner committed
523

524 525
    do
    {
526 527
        ct = 0;
        for (int i = 0; i < m_octants.size(); i++)
528
        {
529
            OctantSharedPtr oct = m_octants[i];
530

531 532
            if (!oct->IsDeltaKnown())
            { // if delta has not been asigned
533
                vector<OctantSharedPtr> known;
534 535
                map<OctantFace, vector<OctantSharedPtr> > nList =
                    oct->GetNeigbours();
536
                map<OctantFace, vector<OctantSharedPtr> >::iterator it;
Michael Turner's avatar
Michael Turner committed
537

538
                for (it = nList.begin(); it != nList.end(); it++)
539
                {
540
                    for (int j = 0; j < it->second.size(); j++)
541
                    {
542
                        if (it->second[j]->IsDeltaKnown())
543 544 545
                        {
                            known.push_back(it->second[j]);
                        }
546
                    }
547
                }
548

549
                if (known.size() > 0)
550 551
                {
                    vector<NekDouble> deltaPrime;
552
                    for (int j = 0; j < known.size(); j++)
553
                    {
554
                        NekDouble r = oct->Distance(known[j]);
555

Michael Turner's avatar
Michael Turner committed
556
                        if (0.199 * r + known[j]->GetDelta() < m_maxDelta)
557
                        {
Michael Turner's avatar
Michael Turner committed
558
                            deltaPrime.push_back(0.199 * r +
559
                                                 known[j]->GetDelta());
560
                        }
561
                        else
562
                        {
563
                            deltaPrime.push_back(m_maxDelta);
564
                        }
565
                    }
566
                    NekDouble min = numeric_limits<double>::max();
567
                    for (int j = 0; j < deltaPrime.size(); j++)
568
                    {
569
                        if (deltaPrime[j] < min)
570
                        {
571
                            min = deltaPrime[j];
572 573
                        }
                    }
574
                    oct->SetDelta(min);
575
                    ct += 1;
576
                    deltaPrime.clear();
577
                }
578
                known.clear();
579
            }
Michael Turner's avatar
Michael Turner committed
580

581 582
            if (oct->GetLocation() == eUnknown)
            { // if the node does not know its location
583
                vector<OctantSharedPtr> known;
584 585
                map<OctantFace, vector<OctantSharedPtr> > nList =
                    oct->GetNeigbours();
586
                map<OctantFace, vector<OctantSharedPtr> >::iterator it;
Michael Turner's avatar
Michael Turner committed
587

588
                for (it = nList.begin(); it != nList.end(); it++)
589
                {
590
                    for (int j = 0; j < it->second.size(); j++)
591
                    {
592
                        if (it->second[j]->GetLocation() != eUnknown)
593 594 595
                        {
                            known.push_back(it->second[j]);
                        }
596
                    }
597
                }
598

599
                if (known.size() > 0)
600
                {
601
                    vector<OctantSharedPtr> isNotOnBound;
602
                    for (int j = 0; j < known.size(); j++)
603
                    {
604
                        if (known[j]->GetLocation() != eOnBoundary)
605
                        {
606
                            isNotOnBound.push_back(known[j]);
607 608
                        }
                    }
609

610
                    if (isNotOnBound.size() > 0)
611
                    {
612
                        oct->SetLocation(isNotOnBound[0]->GetLocation());
613 614 615
                    }
                    else
                    {
616
                        NekDouble dist = numeric_limits<double>::max();
617 618 619

                        OctantSharedPtr closest;

Michael Turner's avatar
Michael Turner committed
620
                        bool f = false;
621
                        for (int j = 0; j < known.size(); j++)
622
                        {
623
                            if (oct->Distance(known[j]) < dist)
624
                            {
625
                                closest = known[j];
626
                                dist    = oct->Distance(known[j]);
mike turner's avatar
mike turner committed
627
                                f       = true;
628 629
                            }
                        }
mike turner's avatar
mike turner committed
630
                        ASSERTL0(f, "closest never set");
631

Michael Turner's avatar
Michael Turner committed
632
                        SPBaseSharedPtr sp = closest->GetABoundPoint();
633

Michael Turner's avatar
Michael Turner committed
634
                        Array<OneD, NekDouble> octloc, sploc, vec(3), uv, N;
635
                        int surf;
Michael Turner's avatar
Michael Turner committed
636
                        sp->GetCAD(surf, uv);
637
                        N = m_mesh->m_cad->GetSurf(surf)->N(uv);
638 639

                        octloc = oct->GetLoc();
Michael Turner's avatar
Michael Turner committed
640
                        sploc  = sp->GetLoc();
641

Michael Turner's avatar
Michael Turner committed
642 643 644
                        vec[0] = octloc[0] - sploc[0];
                        vec[1] = octloc[1] - sploc[1];
                        vec[2] = octloc[2] - sploc[2];
645

646 647
                        NekDouble dot =
                            vec[0] * N[0] + vec[1] * N[1] + vec[2] * N[2];
648

649
                        if (dot <= 0.0)
650
                        {
651
                            oct->SetLocation(eOutside);
652
                            ct += 1;
653 654 655
                        }
                        else
                        {
656
                            oct->SetLocation(eInside);
657
                            ct += 1;
658
                        }
659 660
                    }
                }
661
                known.clear();
662
            }
663
        }
Michael Turner's avatar
Michael Turner committed
664

665
    } while (ct > 0);
Michael Turner's avatar
Michael Turner committed
666

667
    for (int i = 0; i < m_octants.size(); i++)
668
    {
669 670 671 672
        if (!m_octants[i]->IsDeltaKnown())
        {
            m_octants[i]->SetDelta(m_maxDelta);
        }
673
    }
674
}
675

676
void Octree::SmoothAllOctants()
677
{
678 679
    // until no more changes occur smooth the mesh specification between all
    // m_octants not particualrly strictly
680
    int ct = 0;