Octree.cpp 35.1 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/CADSurf.h>
38
#include <NekMeshUtils/CADSystem/CADCurve.h>
39
#include <NekMeshUtils/Module/Module.h>
40

Michael Turner's avatar
Michael Turner committed
41
#include <LibUtilities/BasicUtils/ParseUtils.hpp>
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);

Michael Turner's avatar
Michael Turner committed
78 79

    //begin recersive subdivision
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
    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;
    }
101 102
}

103 104
NekDouble Octree::Query(Array<OneD, NekDouble> loc)
{
105
    // starting at master octant 0 move through succsesive m_octants which
Michael Turner's avatar
Michael Turner committed
106
    // contain the point loc until a leaf is found
Michael Turner's avatar
Michael Turner committed
107 108 109 110 111 112 113 114 115 116
    //first search through sourcepoints

    for(int i = 0; i < m_lsources.size(); i++)
    {
        if(m_lsources[i].withinRange(loc))
        {
            return m_lsources[i].delta;
        }
    }

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

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

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

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

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

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

189
void Octree::WriteOctree(string nm)
190
{
191 192 193 194 195
    MeshSharedPtr oct = boost::shared_ptr<Mesh>(new Mesh());
    oct->m_expDim     = 3;
    oct->m_spaceDim   = 3;
    oct->m_nummode    = 2;

196
    for (int i = 0; i < m_octants.size(); i++)
197
    {
198
        /*if(m_octants[i]->GetLocation() != eOnBoundary)
199 200 201
        {
            continue;
        }*/
202 203 204

        vector<NodeSharedPtr> ns(8);

205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243
        ns[0] = boost::shared_ptr<Node>(new Node(0,
                                                 m_octants[i]->FX(eBack),
                                                 m_octants[i]->FX(eDown),
                                                 m_octants[i]->FX(eRight)));

        ns[1] = boost::shared_ptr<Node>(new Node(0,
                                                 m_octants[i]->FX(eForward),
                                                 m_octants[i]->FX(eDown),
                                                 m_octants[i]->FX(eRight)));

        ns[2] = boost::shared_ptr<Node>(new Node(0,
                                                 m_octants[i]->FX(eForward),
                                                 m_octants[i]->FX(eUp),
                                                 m_octants[i]->FX(eRight)));

        ns[3] = boost::shared_ptr<Node>(new Node(0,
                                                 m_octants[i]->FX(eBack),
                                                 m_octants[i]->FX(eUp),
                                                 m_octants[i]->FX(eRight)));

        ns[4] = boost::shared_ptr<Node>(new Node(0,
                                                 m_octants[i]->FX(eBack),
                                                 m_octants[i]->FX(eDown),
                                                 m_octants[i]->FX(eLeft)));

        ns[5] = boost::shared_ptr<Node>(new Node(0,
                                                 m_octants[i]->FX(eForward),
                                                 m_octants[i]->FX(eDown),
                                                 m_octants[i]->FX(eLeft)));

        ns[6] = boost::shared_ptr<Node>(new Node(0,
                                                 m_octants[i]->FX(eForward),
                                                 m_octants[i]->FX(eUp),
                                                 m_octants[i]->FX(eLeft)));

        ns[7] = boost::shared_ptr<Node>(new Node(0,
                                                 m_octants[i]->FX(eBack),
                                                 m_octants[i]->FX(eUp),
                                                 m_octants[i]->FX(eLeft)));
244 245 246

        vector<int> tags;
        tags.push_back(0);
247
        ElmtConfig conf(LibUtilities::eHexahedron, 1, false, false);
248
        ElementSharedPtr E = GetElementFactory().CreateInstance(
249
            LibUtilities::eHexahedron, conf, ns, tags);
250
        oct->m_element[3].push_back(E);
251
    }
Michael Turner's avatar
Michael Turner committed
252

253 254 255 256 257 258 259 260 261
    ModuleSharedPtr mod = GetModuleFactory().CreateInstance(
        ModuleKey(eOutputModule, "xml"), oct);
    mod->RegisterConfig("outfile", nm);
    mod->ProcessVertices();
    mod->ProcessEdges();
    mod->ProcessFaces();
    mod->ProcessElements();
    mod->ProcessComposites();
    mod->Process();
Michael Turner's avatar
Michael Turner committed
262
}
Michael Turner's avatar
Michael Turner committed
263

264
void Octree::SubDivide()
Michael Turner's avatar
Michael Turner committed
265
{
266
    bool repeat;
267 268
    int ct = 1;

269 270 271
    m_numoct = 1;
    m_masteroct->Subdivide(m_masteroct, m_numoct);

272
    if (m_mesh->m_verbose)
273 274 275
    {
        cout << "\tSubdivide iteration: ";
    }
Michael Turner's avatar
Michael Turner committed
276 277 278

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

292
        VerifyNeigbours();
293

294 295 296 297 298 299 300 301
        // 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.
302 303
        vector<vector<OctantSharedPtr> > dividelist;
        set<int> inlist;
304
        // build initial list
305
        {
306
            vector<OctantSharedPtr> sublist;
307
            for (int i = 0; i < m_octants.size(); i++)
308
            {
309 310
                if (m_octants[i]->NeedDivide() &&
                    m_octants[i]->DX() / 4.0 > m_minDelta)
311
                {
312 313
                    sublist.push_back(m_octants[i]);
                    inlist.insert(m_octants[i]->GetId());
314 315
                    repeat = true; // if there is a subdivision this whole
                                   // process needs to be repeated
316 317
                }
            }
318
            dividelist.push_back(sublist);
319
        }
320
        // then loop over building sublists until no more are required
321
        int ct2 = 0;
322
        while (true)
Michael Turner's avatar
Michael Turner committed
323
        {
324
            ct2++;
325 326 327
            vector<OctantSharedPtr> newsublist,
                previouslist = dividelist.back();
            for (int i = 0; i < previouslist.size(); i++)
Michael Turner's avatar
Michael Turner committed
328
            {
329 330
                map<OctantFace, vector<OctantSharedPtr> > nlist =
                    previouslist[i]->GetNeigbours();
331
                map<OctantFace, vector<OctantSharedPtr> >::iterator it;
332
                for (it = nlist.begin(); it != nlist.end(); it++)
Michael Turner's avatar
Michael Turner committed
333
                {
334
                    for (int j = 0; j < it->second.size(); j++)
Michael Turner's avatar
Michael Turner committed
335
                    {
336
                        if (previouslist[i]->DX() < it->second[j]->DX())
Michael Turner's avatar
Michael Turner committed
337
                        {
338 339 340
                            set<int>::iterator s =
                                inlist.find(it->second[j]->GetId());
                            if (s == inlist.end())
341 342
                            {
                                inlist.insert(it->second[j]->GetId());
343
                                newsublist.push_back(it->second[j]);
344
                            }
Michael Turner's avatar
Michael Turner committed
345 346 347 348
                        }
                    }
                }
            }
349
            if (newsublist.size() == 0)
350 351 352 353 354 355 356
            {
                break;
            }
            else
            {
                dividelist.push_back(newsublist);
            }
Michael Turner's avatar
Michael Turner committed
357
        }
358 359

        vector<vector<OctantSharedPtr> >::reverse_iterator rit;
360
        for (rit = dividelist.rbegin(); rit != dividelist.rend(); rit++)
361
        {
362
            vector<OctantSharedPtr> currentlist = *rit;
363
            for (int i = 0; i < currentlist.size(); i++)
364
            {
365
                currentlist[i]->Subdivide(currentlist[i], m_numoct);
366
            }
367
        }
368
    } while (repeat);
369

370
    if (m_mesh->m_verbose)
371 372 373
    {
        cout << endl;
    }
374
}
375

376
bool Octree::VerifyNeigbours()
377
{
Michael Turner's avatar
Michael Turner committed
378 379 380 381
    // 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
382
    bool error = false;
383
    for (int i = 0; i < m_octants.size(); i++)
384 385
    {
        bool valid = true;
386 387
        map<OctantFace, vector<OctantSharedPtr> > nlist =
            m_octants[i]->GetNeigbours();
388
        map<OctantFace, vector<OctantSharedPtr> >::iterator it;
389
        for (it = nlist.begin(); it != nlist.end(); it++)
390
        {
391
            if (it->second.size() == 0)
392 393 394 395 396
            {
                NekDouble expectedfx;
                switch (it->first)
                {
                    case eUp:
397
                        expectedfx = m_centroid[1] + m_dim;
398 399
                        break;
                    case eDown:
400
                        expectedfx = m_centroid[1] - m_dim;
401 402
                        break;
                    case eLeft:
403
                        expectedfx = m_centroid[2] + m_dim;
404 405
                        break;
                    case eRight:
406
                        expectedfx = m_centroid[2] - m_dim;
407 408
                        break;
                    case eForward:
409
                        expectedfx = m_centroid[0] + m_dim;
410 411
                        break;
                    case eBack:
412
                        expectedfx = m_centroid[0] - m_dim;
413 414
                        break;
                }
415
                if (fabs(m_octants[i]->FX(it->first) - expectedfx) > 1E-6)
416 417 418
                {
                    valid = false;
                    cout << "wall neigbour error" << endl;
419 420
                    cout << expectedfx << " " << m_octants[i]->FX(it->first)
                         << " " << it->first << endl;
421 422
                }
            }
423
            else if (it->second.size() == 1)
424
            {
425 426
                if (!(m_octants[i]->DX() == it->second[0]->DX() ||
                      it->second[0]->DX() == 2.0 * m_octants[i]->DX()))
427 428 429
                {
                    valid = false;
                    cout << " 1 neigbour error" << endl;
430 431
                    cout << m_octants[i]->DX() << " " << it->second[0]->DX()
                         << endl;
432 433
                }
            }
434
            else if (it->second.size() == 4)
435
            {
436
                if (!(m_octants[i]->DX() / 2.0 == it->second[0]->DX()))
437 438 439
                {
                    valid = false;
                    cout << "4 neibour error" << endl;
440 441
                    cout << m_octants[i]->DX() << " " << it->second[0]->DX()
                         << endl;
442
                }
443
            }
444
        }
445
        if (!valid)
446 447 448 449
        {
            error = true;
            cout << "invalid neigbour config" << endl;
        }
450
    }
451
    return !error;
452 453
}

454
void Octree::SmoothSurfaceOctants()
455
{
456 457 458
    // 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
459
    int ct = 0;
Michael Turner's avatar
Michael Turner committed
460

461 462
    do
    {
463 464
        ct = 0;
        for (int i = 0; i < m_octants.size(); i++)
465
        {
466
            OctantSharedPtr oct = m_octants[i];
467

468
            if (oct->IsDeltaKnown())
469
            {
470
                vector<OctantSharedPtr> check;
471
                map<OctantFace, vector<OctantSharedPtr> > nList =
472
                    oct->GetNeigbours();
473
                map<OctantFace, vector<OctantSharedPtr> >::iterator it;
Michael Turner's avatar
Michael Turner committed
474

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

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

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

511 512
void Octree::PropagateDomain()
{
513 514 515 516 517
    // 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
518

519 520
    do
    {
521 522
        ct = 0;
        for (int i = 0; i < m_octants.size(); i++)
523
        {
524
            OctantSharedPtr oct = m_octants[i];
525

526 527
            if (!oct->IsDeltaKnown())
            { // if delta has not been asigned
528
                vector<OctantSharedPtr> known;
529 530
                map<OctantFace, vector<OctantSharedPtr> > nList =
                    oct->GetNeigbours();
531
                map<OctantFace, vector<OctantSharedPtr> >::iterator it;
Michael Turner's avatar
Michael Turner committed
532

533
                for (it = nList.begin(); it != nList.end(); it++)
534
                {
535
                    for (int j = 0; j < it->second.size(); j++)
536
                    {
537
                        if (it->second[j]->IsDeltaKnown())
538 539 540
                        {
                            known.push_back(it->second[j]);
                        }
541
                    }
542
                }
543

544
                if (known.size() > 0)
545 546
                {
                    vector<NekDouble> deltaPrime;
547
                    for (int j = 0; j < known.size(); j++)
548
                    {
549
                        NekDouble r = oct->Distance(known[j]);
550

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

576 577
            if (oct->GetLocation() == eUnknown)
            { // if the node does not know its location
578
                vector<OctantSharedPtr> known;
579 580
                map<OctantFace, vector<OctantSharedPtr> > nList =
                    oct->GetNeigbours();
581
                map<OctantFace, vector<OctantSharedPtr> >::iterator it;
Michael Turner's avatar
Michael Turner committed
582

583
                for (it = nList.begin(); it != nList.end(); it++)
584
                {
585
                    for (int j = 0; j < it->second.size(); j++)
586
                    {
587
                        if (it->second[j]->GetLocation() != eUnknown)
588 589 590
                        {
                            known.push_back(it->second[j]);
                        }
591
                    }
592
                }
593

594
                if (known.size() > 0)
595
                {
596
                    vector<OctantSharedPtr> isNotOnBound;
597
                    for (int j = 0; j < known.size(); j++)
598
                    {
599
                        if (known[j]->GetLocation() != eOnBoundary)
600
                        {
601
                            isNotOnBound.push_back(known[j]);
602 603
                        }
                    }
604

605
                    if (isNotOnBound.size() > 0)
606
                    {
607
                        oct->SetLocation(isNotOnBound[0]->GetLocation());
608 609 610
                    }
                    else
                    {
611
                        NekDouble dist = numeric_limits<double>::max();
612 613 614

                        OctantSharedPtr closest;

Michael Turner's avatar
Michael Turner committed
615
                        bool f = false;
616
                        for (int j = 0; j < known.size(); j++)
617
                        {
618
                            if (oct->Distance(known[j]) < dist)
619
                            {
620
                                closest = known[j];
621
                                dist    = oct->Distance(known[j]);
Michael Turner's avatar
Michael Turner committed
622
                                f = true;
623 624
                            }
                        }
Michael Turner's avatar
Michael Turner committed
625
                        ASSERTL0(f,"closest never set");
626

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

Michael Turner's avatar
Michael Turner committed
629
                        Array<OneD, NekDouble> octloc, sploc, vec(3), uv, N;
630
                        int surf;
Michael Turner's avatar
Michael Turner committed
631
                        sp->GetCAD(surf, uv);
632
                        N = m_mesh->m_cad->GetSurf(surf)->N(uv);
633 634

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

Michael Turner's avatar
Michael Turner committed
637 638 639
                        vec[0] = octloc[0] - sploc[0];
                        vec[1] = octloc[1] - sploc[1];
                        vec[2] = octloc[2] - sploc[2];
640

641 642
                        NekDouble dot =
                            vec[0] * N[0] + vec[1] * N[1] + vec[2] * N[2];
643

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

660
    } while (ct > 0);
Michael Turner's avatar
Michael Turner committed
661

662
    for (int i = 0; i < m_octants.size(); i++)
663
    {
664 665
        ASSERTL0(m_octants[i]->IsDeltaKnown(),
                 "does not know delta after propergation");
666
    }
667
}
668

669
void Octree::SmoothAllOctants()
670
{
671 672
    // until no more changes occur smooth the mesh specification between all
    // m_octants not particualrly strictly
673
    int ct = 0;
Michael Turner's avatar
Michael Turner committed
674

675 676
    do
    {
677 678
        ct = 0;
        for (int i = 0; i < m_octants.size(); i++)
679
        {
680
            OctantSharedPtr oct = m_octants[i];
681

682
            vector<OctantSharedPtr> check;
683 684
            map<OctantFace, vector<OctantSharedPtr> > nList =
                oct->GetNeigbours();
685
            map<OctantFace, vector<OctantSharedPtr> >::iterator it;
686
            for (it = nList.begin(); it != nList.end(); it++)
687
            {
688
                for (int j = 0; j < it->second.size(); j++)
689
                {
690 691
                    if (it->second[j]->GetDelta() < oct->GetDelta() &&
                        ddx(oct, it->second[j]) > 0.2)
692
                    {
693
                        check.push_back(it->second[j]);
694
                    }
695
                }
696
            }
Michael Turner's avatar
Michael Turner committed
697

698
            if (check.size() > 0)
699 700
            {
                NekDouble deltaSM = numeric_limits<double>::max();
701
                for (int j = 0; j < check.size(); j++)
702
                {
703
                    NekDouble r = oct->Distance(check[j]);
Michael Turner's avatar
Michael Turner committed
704

705
                    if (0.199 * r + check[j]->GetDelta() < deltaSM)
706 707
                    {
                        deltaSM = 0.199 * r + check[j]->GetDelta();
708 709
                    }
                }
710
                oct->SetDelta(deltaSM);
711
                ct += 1;
712
            }
713
        }
Michael Turner's avatar
Michael Turner committed
714

715
    } while (ct > 0);
716 717
}

718
int Octree::CountElemt()
719
{
720 721
    // by considering the volume of a tet evaluate the number of elements in the
    // mesh
Michael Turner's avatar
Michael Turner committed
722

723
    NekDouble total = 0.0;
Michael Turner's avatar
Michael Turner committed
724

725
    Array<OneD, NekDouble> boundingBox = m_mesh->m_cad->GetBoundingBox();
726

727
    for (int i = 0; i < m_octants.size(); i++)
728
    {
729
        OctantSharedPtr oct = m_octants[i];
730
        if (oct->GetLocation() == eInside)
731
        {
732 733 734
            total += 8.0 * oct->DX() * oct->DX() * oct->DX() /
                     (oct->GetDelta() * oct->GetDelta() * oct->GetDelta() /
                      6.0 / sqrt(2));
735
        }
736
        else if (oct->GetLocation() == eOnBoundary)
737
        {
738
            NekDouble vol = 1.0;
739 740
            if (oct->FX(eBack) < boundingBox[1] &&
                oct->FX(eForward) > boundingBox[0])
741
            {
742
                // then there is some over lap in x
743
                NekDouble min, max;
744
                if (oct->FX(eBack) > boundingBox[0])
745 746 747 748 749 750 751
                {
                    min = oct->FX(eBack);
                }
                else
                {
                    min = boundingBox[0];
                }
752
                if (boundingBox[1] < oct->FX(eForward))
753 754 755 756 757 758 759 760 761 762 763 764 765 766
                {
                    max = boundingBox[1];
                }
                else
                {
                    max = oct->FX(eForward);
                }
                vol *= (max - min);
            }
            else
            {
                vol *= 0.0;
            }

767 768
            if (oct->FX(eDown) < boundingBox[3] &&
                oct->FX(eUp) > boundingBox[2])
769
            {
770
                // then there is some over lap in x
771
                NekDouble min, max;
772
                if (oct->FX(eDown) > boundingBox[2])
773 774 775 776 777 778 779
                {
                    min = oct->FX(eDown);
                }
                else
                {
                    min = boundingBox[2];
                }
780
                if (boundingBox[3] < oct->FX(eUp))
781 782 783 784 785 786 787 788 789 790 791 792 793
                {
                    max = boundingBox[3];
                }
                else
                {
                    max = oct->FX(eUp);
                }
                vol *= (max - min);
            }
            else
            {
                vol *= 0.0;
            }
794

795 796
            if (oct->FX(eRight) < boundingBox[5] &&
                oct->FX(eLeft) > boundingBox[4])
797
            {
798
                // then there is some over lap in x
799
                NekDouble min, max;
800
                if (oct->FX(eRight) > boundingBox[4])
801 802 803 804 805 806 807
                {
                    min = oct->FX(eRight);
                }
                else
                {
                    min = boundingBox[4];
                }
808
                if (boundingBox[5] < oct->FX(eLeft))
809 810 811 812 813 814 815 816 817 818 819 820 821
                {
                    max = boundingBox[5];
                }
                else
                {
                    max = oct->FX(eLeft);
                }
                vol *= (max - min);
            }
            else
            {
                vol *= 0.0;
            }
822 823
            total += vol / 2.0 / (oct->GetDelta() * oct->GetDelta() *
                                  oct->GetDelta() / 6.0 / sqrt(2));
824
        }
825 826
    }

827
    return int(total);
828
}
829

830
void Octree::CompileSourcePointList()
831
{
832 833 834 835 836 837

    for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); i++)
    {
        CADCurveSharedPtr curve = m_mesh->m_cad->GetCurve(i);
        Array<OneD, NekDouble> bds = curve->Bounds();
        int samples = 100;
Michael Turner's avatar
Michael Turner committed
838
        NekDouble dt      = (bds[1] - bds[0]) / (samples + 1);
839
        for (int j = 1; j < samples -1; j++) //dont want first and last point
840
        {
Michael Turner's avatar