ProcessBL.cpp 38.4 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
////////////////////////////////////////////////////////////////////////////////
//
//  File: ProcessBL.cpp
//
//  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.
//
32
//  Description: Refine prismatic or quadrilateral boundary layer elements.
33 34 35
//
////////////////////////////////////////////////////////////////////////////////

36 37
#include <string>

38
#include <LibUtilities/BasicUtils/ParseUtils.h>
39 40 41
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Foundations/BLPoints.h>
#include <LibUtilities/Foundations/ManagerAccess.h>
42
#include <LibUtilities/Interpreter/AnalyticExpressionEvaluator.hpp>
43
#include <LocalRegions/HexExp.h>
44
#include <LocalRegions/PrismExp.h>
45
#include <LocalRegions/QuadExp.h>
46

47
#include <NekMeshUtils/MeshElements/Element.h>
48 49 50
#include "ProcessBL.h"

using namespace std;
51
using namespace Nektar::NekMeshUtils;
52

53 54
namespace Nektar
{
55 56
namespace Utilities
{
57 58 59
ModuleKey ProcessBL::className = GetModuleFactory().RegisterCreatorFunction(
    ModuleKey(eProcessModule, "bl"),
    ProcessBL::create,
60
    "Refines a prismatic or quadrilateral boundary layer.");
61 62 63

int **helper2d(int lda, int arr[][2])
{
64
    int **ret = new int *[lda];
65
    for (int i = 0; i < lda; ++i)
66
    {
67
        ret[i]    = new int[2];
68 69 70 71 72
        ret[i][0] = arr[i][0];
        ret[i][1] = arr[i][1];
    }
    return ret;
}
73

74 75
int **helper2d(int lda, int arr[][4])
{
76
    int **ret = new int *[lda];
77 78
    for (int i = 0; i < lda; ++i)
    {
79
        ret[i]    = new int[4];
80 81 82 83 84 85 86
        ret[i][0] = arr[i][0];
        ret[i][1] = arr[i][1];
        ret[i][2] = arr[i][2];
        ret[i][3] = arr[i][3];
    }
    return ret;
}
87

88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
struct SplitMapHelper
{
    int size;
    int layerOff;
    int *edge;
    int *offset;
    int *inc;
    int **conn;
    int bfacesSize;
    int *bfaces;
};

struct SplitEdgeHelper
{
    int size;
    int *edge;
    int **edgeVert;
    int *offset;
    int *inc;
};

ProcessBL::ProcessBL(MeshSharedPtr m) : ProcessModule(m)
{
    // BL mesh configuration.
112 113 114 115 116 117 118 119
    m_config["layers"] =
        ConfigOption(false, "2", "Number of layers to refine.");
    m_config["nq"] =
        ConfigOption(false, "5", "Number of points in high order elements.");
    m_config["surf"] =
        ConfigOption(false, "", "Tag identifying surface connected to prism.");
    m_config["r"] =
        ConfigOption(false, "2.0", "Ratio to use in geometry progression.");
120
}
121

122 123 124
ProcessBL::~ProcessBL()
{
}
125

126 127 128 129
void ProcessBL::Process()
{
    if (m_mesh->m_verbose)
    {
130
        cout << "ProcessBL: Refining boundary layer..." << endl;
131
    }
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
    int dim = m_mesh->m_expDim;
    switch (dim)
    {
        case 2:
        {
            BoundaryLayer2D();
        }
        break;

        case 3:
        {
            BoundaryLayer3D();
        }
        break;

        default:
            ASSERTL0(0,"Dimension not supported");
        break;
    }

    // Re-process mesh to eliminate duplicate vertices and edges.
    ProcessVertices();
    ProcessEdges();
    ProcessFaces();
    ProcessElements();
    ProcessComposites();
}

void ProcessBL::BoundaryLayer2D()
{
    int nodeId  = m_mesh->m_vertexSet.size();
    int nl      = m_config["layers"].as<int>();
    int nq      = m_config["nq"].    as<int>();

Michael Turner's avatar
mac  
Michael Turner committed
166
    // determine if geometric ratio is string or a constant.
167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
    LibUtilities::AnalyticExpressionEvaluator rEval;
    NekDouble r             =  1;
    int       rExprId       = -1;
    bool      ratioIsString = false;

    if (m_config["r"].isType<NekDouble>())
    {
        r = m_config["r"].as<NekDouble>();
    }
    else
    {
        std::string rstr = m_config["r"].as<string>();
        rExprId = rEval.DefineFunction("x y z", rstr);
        ratioIsString = true;
    }

    // Default PointsType.
    LibUtilities::PointsType pt = LibUtilities::eGaussLobattoLegendre;

    // Map which takes element ID to edge on surface. This enables
    // splitting to occur in either y-direction of the prism.
    map<int, int> splitEls;

    // edgeMap associates geometry edge IDs to the (nl+1) vertices which
    // are generated along that edge when a prism is split, and is used
    // to avoid generation of duplicate vertices. It is stored as an
    // unordered map for speed.
194
    std::unordered_map<int, vector<NodeSharedPtr> > edgeMap;
195 196 197 198 199

    string surf = m_config["surf"].as<string>();
    if (surf.size() > 0)
    {
        vector<unsigned int> surfs;
200
        ParseUtils::GenerateSeqVector(surf, surfs);
201 202 203 204 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 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290
        sort(surfs.begin(), surfs.end());

        // If surface is defined, process list of elements to find those
        // that are connected to it.
        for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
        {
            ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
            int nSurf = el->GetEdgeCount();

            for (int j = 0; j < nSurf; ++j)
            {
                int bl = el->GetBoundaryLink(j);
                if (bl == -1)
                {
                    continue;
                }

                ElementSharedPtr bEl  = m_mesh->m_element[m_mesh->m_expDim-1][bl];
                vector<int>      tags = bEl->GetTagList();
                vector<int>      inter;

                sort(tags.begin(), tags.end());
                set_intersection(surfs.begin(), surfs.end(),
                                 tags .begin(), tags .end(),
                                 back_inserter(inter));
                ASSERTL0(inter.size() <= 1,
                         "Intersection of surfaces wrong");

                if (inter.size() == 1)
                {
                    if (el->GetConf().m_e != LibUtilities::eQuadrilateral)
                    {
                        cerr << "WARNING: Found non-quad element "
                             << "to split in surface " << surf
                             << "; ignoring" << endl;
                        continue;
                    }

                    if (splitEls.count(el->GetId()) > 0)
                    {
                        cerr << "WARNING: quad already found; "
                             << "ignoring" << endl;
                        continue;
                    }

                    splitEls[el->GetId()] = j;
                }
            }
        }
    }
    else
    {
        ASSERTL0(false, "Surface must be specified.");
    }

    if (splitEls.size() == 0)
    {
        cerr << "WARNING: No elements detected to split." << endl;
        return;
    }

    // Erase all elements from the element list. Elements will be
    // re-added as they are split.
    vector<ElementSharedPtr> el = m_mesh->m_element[m_mesh->m_expDim];
    m_mesh->m_element[m_mesh->m_expDim].clear();

    // Iterate over list of elements of expansion dimension.
    for (int i = 0; i < el.size(); ++i)
    {


        if (splitEls.count(el[i]->GetId()) == 0)
        {
            m_mesh->m_element[m_mesh->m_expDim].push_back(el[i]);
            continue;
        }

        // Find other boundary faces if any
        std::map<int, int> bLink;
        for (int j = 0; j < 4; j++)
        {
            int bl = el[i]->GetBoundaryLink(j);
            if ( (bl != -1) && ( j != splitEls[el[i]->GetId()]) )
            {
                bLink[j] = bl;
            }
        }

        // Get elemental geometry object.
        SpatialDomains::QuadGeomSharedPtr geom =
291
            std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
292 293
                el[i]->GetGeom(m_mesh->m_spaceDim));

Michael Turner's avatar
mac  
Michael Turner committed
294
        // Determine whether to use reverse points.
295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341
        // (if edges 1 or 2 are on the surface)
        LibUtilities::PointsType t = ( (splitEls[el[i]->GetId()]+1) %4) < 2 ?
            LibUtilities::eBoundaryLayerPoints :
            LibUtilities::eBoundaryLayerPointsRev;


        if(ratioIsString) // determine value of r base on geom
        {
            NekDouble x,y,z;
            NekDouble x1,y1,z1;
            int nverts = geom->GetNumVerts();

            x = y = z = 0.0;

            for(int i = 0; i < nverts; ++i)
            {
                geom->GetVertex(i)->GetCoords(x1,y1,z1);
                x += x1; y += y1; z += z1;
            }
            x /= (NekDouble) nverts;
            y /= (NekDouble) nverts;
            z /= (NekDouble) nverts;
            r = rEval.Evaluate(rExprId,x,y,z,0.0);
        }

        // Create basis.
        LibUtilities::BasisKey B0(
            LibUtilities::eModified_A, nq,
            LibUtilities::PointsKey(nq,pt));
        LibUtilities::BasisKey B1(
            LibUtilities::eModified_A, 2,
            LibUtilities::PointsKey(nl+1, t, r));

        // Create local region.
        LocalRegions::QuadExpSharedPtr q;
        if (splitEls[el[i]->GetId()] % 2)
        {
            q = MemoryManager<LocalRegions::QuadExp>::AllocateSharedPtr(
                B1,B0,geom);
        }
        else
        {
            q = MemoryManager<LocalRegions::QuadExp>::AllocateSharedPtr(
                B0,B1,geom);
        }

        // Grab co-ordinates.
Julian Marcon's avatar
Julian Marcon committed
342 343
        Array<OneD, NekDouble> x(nq*(nl+1), 0.0);
        Array<OneD, NekDouble> y(nq*(nl+1), 0.0);
Julian Marcon's avatar
Julian Marcon committed
344
        Array<OneD, NekDouble> z(nq*(nl+1), 0.0);
345 346 347 348 349 350 351 352 353 354 355 356
        q->GetCoords(x,y,z);

        vector<vector<NodeSharedPtr> > edgeNodes(2);

        // Loop over edges to be split.
        for (int j = 0; j < 2; ++j)
        {
            int locEdge = (splitEls[el[i]->GetId()]+1+2*j)%4;
            int edgeId  = el[i]->GetEdge(locEdge)->m_id;

            // Determine whether we have already generated vertices
            // along this edge.
357
            auto eIt = edgeMap.find(edgeId);
358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373

            if (eIt == edgeMap.end())
            {
                // If not then resize storage to hold new points.
                edgeNodes[j].resize(nl+1);

                // Re-use existing vertices at endpoints of edge to
                // avoid duplicating the existing vertices.
                edgeNodes[j][0]  = el[i]->GetVertex(locEdge);
                edgeNodes[j][nl] = el[i]->GetVertex((locEdge+1)%4);

                 // Variable geometric ratio
                if(ratioIsString)
                {
                    NekDouble x0,y0;
                    NekDouble x1,y1;
Michael Turner's avatar
mac  
Michael Turner committed
374
                    NekDouble xm,ym,zm=0.0;
375 376 377 378 379 380 381 382 383 384 385 386 387 388 389

                    // -> Find edge end and mid points
                    x0 = edgeNodes[j][0]->m_x;
                    y0 = edgeNodes[j][0]->m_y;

                    x1 = edgeNodes[j][nl]->m_x;
                    y1 = edgeNodes[j][nl]->m_y;

                    xm = 0.5*(x0+x1);
                    ym = 0.5*(y0+y1);

                    // evaluate r factor based on mid point value
                    NekDouble rnew;
                    rnew = rEval.Evaluate(rExprId,xm,ym,zm,0.0);

Michael Turner's avatar
mac  
Michael Turner committed
390
                    // Get basis with new r;
391 392 393 394 395 396 397
                    t =  (j==0) ? LibUtilities::eBoundaryLayerPoints :
                                 LibUtilities::eBoundaryLayerPointsRev;
                    LibUtilities::PointsKey Pkey(nl+1, t, rnew);
                    LibUtilities::PointsSharedPtr newP
                        = LibUtilities::PointsManager()[Pkey];

                    const Array<OneD, const NekDouble> z = newP->GetZ();
398

399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428
                    // Create new interior nodes based on this new blend
                    for (int k = 1; k < nl; ++k)
                    {
                        xm = 0.5*(1+z[k])*(x1-x0) + x0;
                        ym = 0.5*(1+z[k])*(y1-y0) + y0;
                        zm = 0.0;
                        edgeNodes[j][k] = NodeSharedPtr(
                                new Node(nodeId++, xm,ym,zm));
                    }
                }
                else
                {
                    // Create new interior nodes.
                    int pos;
                    for (int k = 1; k < nl; ++k)
                    {
                        switch (locEdge)
                        {
                            case 0:
                                pos = k;
                                break;
                            case 1:
                                pos = nq -1 + k*nq;
                                break;
                            case 2:
                                pos = nq*(nl+1) -1 - k;
                                break;
                            case 3:
                                pos = nq*nl - k*nq;
                                break;
429 430 431
                            default:
                                ASSERTL0(0,"Quad edge should be < 4.");
                            break;
432 433 434 435 436 437 438 439 440 441 442 443
                        }
                        edgeNodes[j][k] = NodeSharedPtr(
                            new Node(nodeId++, x[pos], y[pos], z[pos]));
                    }
                }

                // Store these edges in edgeMap.
                edgeMap[edgeId] = edgeNodes[j];
            }
            else
            {
                // Check orientation
444
                if (eIt->second[0] == el[i]->GetVertex(locEdge))
445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530
                {
                    // Same orientation: copy nodes
                    edgeNodes[j] = eIt->second;
                }
                else
                {
                    // Reversed orientation: copy in reversed order
                    edgeNodes[j].resize(nl+1);
                    for (int k = 0; k < nl+1; ++k)
                    {
                        edgeNodes[j][k] = eIt->second[nl-k];
                    }
                }

            }
        }

        // Create element layers.
        for (int j = 0; j < nl; ++j)
        {
            // Get corner vertices.
            vector<NodeSharedPtr> nodeList(4);
            switch (splitEls[el[i]->GetId()])
            {
                case 0:
                {
                    nodeList[0] = edgeNodes[1][nl-j  ];
                    nodeList[1] = edgeNodes[0][j  ];
                    nodeList[2] = edgeNodes[0][j+1];
                    nodeList[3] = edgeNodes[1][nl-j-1];
                    break;
                }
                case 1:
                {
                    nodeList[0] = edgeNodes[1][j];
                    nodeList[1] = edgeNodes[1][j+1];
                    nodeList[2] = edgeNodes[0][nl-j-1];
                    nodeList[3] = edgeNodes[0][nl-j];
                    break;
                }
                case 2:
                {
                    nodeList[0] = edgeNodes[0][nl-j  ];
                    nodeList[1] = edgeNodes[1][j  ];
                    nodeList[2] = edgeNodes[1][j+1];
                    nodeList[3] = edgeNodes[0][nl-j-1];
                    break;
                }
                case 3:
                {
                    nodeList[0] = edgeNodes[0][j];
                    nodeList[1] = edgeNodes[0][j+1];
                    nodeList[2] = edgeNodes[1][nl-j-1];
                    nodeList[3] = edgeNodes[1][nl-j];
                    break;
                }
            }
            // Create the element.
            ElmtConfig conf(LibUtilities::eQuadrilateral, 1, true, false, true);
            ElementSharedPtr elmt = GetElementFactory().
                CreateInstance(
              LibUtilities::eQuadrilateral,conf,nodeList,el[i]->GetTagList());

            // Add high order nodes to split edges.
            for (int l = 0; l < 2; ++l)
            {
                int locEdge = (splitEls[el[i]->GetId()]+2*l)%4;
                EdgeSharedPtr HOedge = elmt->GetEdge(
                    locEdge);
                int pos;
                for (int k = 1; k < nq-1; ++k)
                {
                    switch (locEdge)
                    {
                        case 0:
                            pos = j*nq+ k;
                            break;
                        case 1:
                            pos = j+1 + k*(nl+1);
                            break;
                        case 2:
                            pos = (j+1)*nq + (nq-1) - k;
                            break;
                        case 3:
                            pos = (nl+1)*(nq-1) + j - k*(nl+1);
                            break;
531 532 533
                        default:
                            ASSERTL0(0,"Quad edge should be < 4.");
                            break;
534 535 536 537 538 539 540 541
                    }
                    HOedge->m_edgeNodes.push_back(
                        NodeSharedPtr(
                            new Node(nodeId++,x[pos],y[pos],0.0)));
                }
                HOedge->m_curveType = pt;
            }

Michael Turner's avatar
mac  
Michael Turner committed
542
            // Change the elements on the boundary
543
            // to match the layers
544
            for (auto &it : bLink)
545
            {
546 547
                int eid = it.first;
                int bl  = it.second;
548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583

                if (j == 0)
                {
                    // For first layer reuse existing 2D element.
                    ElementSharedPtr e = m_mesh->m_element[m_mesh->m_expDim-1][bl];
                    for (int k = 0; k < 2; ++k)
                    {
                        e->SetVertex(
                            k, nodeList[(eid+k)%4]);
                    }
                }
                else
                {
                    // For all other layers create new element.
                    vector<NodeSharedPtr> qNodeList(2);
                    for (int k = 0; k < 2; ++k)
                    {
                        qNodeList[k] = nodeList[(eid+k)%4];
                    }
                    vector<int> tagBE;
                    tagBE = m_mesh->m_element[m_mesh->m_expDim-1][bl]->GetTagList();
                    ElmtConfig bconf(LibUtilities::eSegment,1,true,true,false);
                    ElementSharedPtr boundaryElmt = GetElementFactory().
                        CreateInstance(LibUtilities::eSegment,bconf,
                                       qNodeList,tagBE);
                    m_mesh->m_element[m_mesh->m_expDim-1].push_back(boundaryElmt);
                }
            }

            m_mesh->m_element[m_mesh->m_expDim].push_back(elmt);
        }
    }
}

void ProcessBL::BoundaryLayer3D()
{
584 585 586 587
    // A set containing all element types which are valid.
    set<LibUtilities::ShapeType> validElTypes;
    validElTypes.insert(LibUtilities::ePrism);
    validElTypes.insert(LibUtilities::eHexahedron);
588

589 590 591
    int nodeId = m_mesh->m_vertexSet.size();
    int nl     = m_config["layers"].as<int>();
    int nq     = m_config["nq"].as<int>();
Dave Moxey's avatar
Dave Moxey committed
592

593 594
    // determine if geometric ratio is string or a constant.
    LibUtilities::AnalyticExpressionEvaluator rEval;
595 596 597
    NekDouble r        = 1;
    int rExprId        = -1;
    bool ratioIsString = false;
Dave Moxey's avatar
Dave Moxey committed
598

599 600 601 602 603 604 605
    if (m_config["r"].isType<NekDouble>())
    {
        r = m_config["r"].as<NekDouble>();
    }
    else
    {
        std::string rstr = m_config["r"].as<string>();
606 607
        rExprId          = rEval.DefineFunction("x y z", rstr);
        ratioIsString    = true;
608
    }
609

610 611
    // Prismatic node -> face map.
    int prismFaceNodes[5][4] = {
612 613 614 615 616 617 618
        {0, 1, 2, 3}, {0, 1, 4, -1}, {1, 2, 5, 4}, {3, 2, 5, -1}, {0, 3, 5, 4}};
    int hexFaceNodes[6][4] = {{0, 1, 2, 3},
                              {0, 1, 5, 4},
                              {1, 2, 6, 5},
                              {3, 2, 6, 7},
                              {0, 3, 7, 4},
                              {4, 5, 6, 7}};
619 620 621 622 623 624 625 626 627
    map<LibUtilities::ShapeType, int **> faceNodeMap;
    faceNodeMap[LibUtilities::ePrism]      = helper2d(5, prismFaceNodes);
    faceNodeMap[LibUtilities::eHexahedron] = helper2d(6, hexFaceNodes);

    // Default PointsType.
    LibUtilities::PointsType pt = LibUtilities::eGaussLobattoLegendre;

    // Map which takes element ID to face on surface. This enables
    // splitting to occur in either y-direction of the prism.
628
    std::unordered_map<int, int> splitEls;
629 630 631 632 633 634 635

    // Set up maps which takes an edge (in nektar++ ordering) and return
    // their offset and stride in the 3d array of collapsed quadrature
    // points. Note that this map includes only the edges that are on
    // the triangular faces as the edges in the normal direction are
    // linear.
    map<LibUtilities::ShapeType, map<int, SplitMapHelper> > splitMap;
636
    int po = nq * (nl + 1);
637 638

    SplitMapHelper splitPrism;
639 640 641 642 643 644 645 646 647 648 649 650 651 652
    int splitMapEdgePrism[6]    = {0, 2, 4, 5, 6, 7};
    int splitMapOffsetPrism[6]  = {0, nq, 0, nq - 1, nq + nq - 1, nq};
    int splitMapIncPrism[6]     = {1, 1, po, po, po, po};
    int splitMapBFacesPrism[3]  = {0, 2, 4};
    int splitMapConnPrism[6][2] = {
        {0, 0}, {1, 0}, {1, 1}, {0, 1}, {2, 0}, {2, 1}};
    splitPrism.size                   = 6;
    splitPrism.layerOff               = nq;
    splitPrism.edge                   = splitMapEdgePrism;
    splitPrism.offset                 = splitMapOffsetPrism;
    splitPrism.inc                    = splitMapIncPrism;
    splitPrism.conn                   = helper2d(6, splitMapConnPrism);
    splitPrism.bfacesSize             = 3;
    splitPrism.bfaces                 = splitMapBFacesPrism;
653 654 655
    splitMap[LibUtilities::ePrism][1] = splitPrism;
    splitMap[LibUtilities::ePrism][3] = splitPrism;

656 657
    int ho = nq * (nq - 1);
    int tl = nq * nq;
658
    SplitMapHelper splitHex0;
659 660 661 662 663 664 665 666 667 668 669 670 671 672 673
    int splitMapEdgeHex0[8]   = {0, 1, 2, 3, 8, 9, 10, 11};
    int splitMapOffsetHex0[8] = {
        0, nq - 1, tl - 1, ho, tl, tl + nq - 1, 2 * tl - 1, tl + ho};
    int splitMapIncHex0[8]     = {1, nq, -1, -nq, 1, nq, -1, -nq};
    int splitMapBFacesHex0[4]  = {1, 2, 3, 4};
    int splitMapConnHex0[8][2] = {
        {0, 0}, {1, 0}, {2, 0}, {3, 0}, {0, 1}, {1, 1}, {2, 1}, {3, 1}};
    splitHex0.size                         = 8;
    splitHex0.layerOff                     = nq * nq;
    splitHex0.edge                         = splitMapEdgeHex0;
    splitHex0.offset                       = splitMapOffsetHex0;
    splitHex0.inc                          = splitMapIncHex0;
    splitHex0.conn                         = helper2d(8, splitMapConnHex0);
    splitHex0.bfacesSize                   = 4;
    splitHex0.bfaces                       = splitMapBFacesHex0;
674 675 676 677 678 679 680 681 682 683 684 685 686 687 688
    splitMap[LibUtilities::eHexahedron][0] = splitHex0;
    splitMap[LibUtilities::eHexahedron][5] = splitHex0;

    // splitEdge enumerates the edges in the standard prism along which
    // new nodes should be generated. These edges are the three between
    // the two triangular faces.
    //
    // edgeVertMap specifies the vertices which comprise those edges in
    // splitEdge; for example splitEdge[0] = 3 which connects vertices 0
    // and 3.
    //
    // edgeOffset holds the offset of each of edges 3, 1 and 8
    // respectively inside the collapsed coordinate system.
    map<LibUtilities::ShapeType, map<int, SplitEdgeHelper> > splitEdge;

689 690 691 692
    int splitPrismEdges[3]       = {3, 1, 8};
    int splitPrismEdgeVert[3][2] = {{0, 3}, {1, 2}, {4, 5}};
    int splitPrismOffset[3]      = {0, nq - 1, nq * (nl + 1) * (nq - 1)};
    int splitPrismInc[3]         = {nq, nq, nq};
693
    SplitEdgeHelper splitPrismEdge;
694 695 696 697 698
    splitPrismEdge.size                = 3;
    splitPrismEdge.edge                = splitPrismEdges;
    splitPrismEdge.edgeVert            = helper2d(3, splitPrismEdgeVert);
    splitPrismEdge.offset              = splitPrismOffset;
    splitPrismEdge.inc                 = splitPrismInc;
699 700 701
    splitEdge[LibUtilities::ePrism][1] = splitPrismEdge;
    splitEdge[LibUtilities::ePrism][3] = splitPrismEdge;

702 703 704 705
    int splitHex0Edges[4]       = {4, 5, 6, 7};
    int splitHex0EdgeVert[4][2] = {{0, 4}, {1, 5}, {2, 6}, {3, 7}};
    int splitHex0Offset[4]      = {0, nq - 1, nq * nq - 1, nq * (nq - 1)};
    int splitHex0Inc[4]         = {nq * nq, nq * nq, nq * nq, nq * nq};
706
    SplitEdgeHelper splitHex0Edge;
707 708 709 710 711
    splitHex0Edge.size                      = 4;
    splitHex0Edge.edge                      = splitHex0Edges;
    splitHex0Edge.edgeVert                  = helper2d(4, splitHex0EdgeVert);
    splitHex0Edge.offset                    = splitHex0Offset;
    splitHex0Edge.inc                       = splitHex0Inc;
712 713 714 715 716 717 718 719 720 721 722 723 724 725
    splitEdge[LibUtilities::eHexahedron][0] = splitHex0Edge;
    splitEdge[LibUtilities::eHexahedron][5] = splitHex0Edge;

    map<LibUtilities::ShapeType, map<int, bool> > revPoints;
    revPoints[LibUtilities::ePrism][1] = true;
    revPoints[LibUtilities::ePrism][3] = false;

    revPoints[LibUtilities::eHexahedron][0] = true;
    revPoints[LibUtilities::eHexahedron][5] = false;

    // edgeMap associates geometry edge IDs to the (nl+1) vertices which
    // are generated along that edge when a prism is split, and is used
    // to avoid generation of duplicate vertices. It is stored as an
    // unordered map for speed.
726
    std::unordered_map<int, vector<NodeSharedPtr> > edgeMap;
727 728 729 730 731

    string surf = m_config["surf"].as<string>();
    if (surf.size() > 0)
    {
        vector<unsigned int> surfs;
732
        ParseUtils::GenerateSeqVector(surf, surfs);
733
        sort(surfs.begin(), surfs.end());
734

735 736 737 738 739
        // If surface is defined, process list of elements to find those
        // that are connected to it.
        for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
        {
            ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
740
            int nSurf           = el->GetFaceCount();
741

742 743 744 745
            for (int j = 0; j < nSurf; ++j)
            {
                int bl = el->GetBoundaryLink(j);
                if (bl == -1)
746
                {
747 748 749
                    continue;
                }

750 751 752 753
                ElementSharedPtr bEl =
                    m_mesh->m_element[m_mesh->m_expDim - 1][bl];
                vector<int> tags = bEl->GetTagList();
                vector<int> inter;
754

755
                sort(tags.begin(), tags.end());
756 757 758 759
                set_intersection(surfs.begin(),
                                 surfs.end(),
                                 tags.begin(),
                                 tags.end(),
760
                                 back_inserter(inter));
761
                ASSERTL0(inter.size() <= 1, "Intersection of surfaces wrong");
762 763 764 765

                if (inter.size() == 1)
                {
                    if (el->GetConf().m_e == LibUtilities::ePrism)
766
                    {
767
                        if (j % 2 == 0)
768
                        {
769 770 771
                            cerr << "WARNING: Found quadrilateral face " << j
                                 << " on surface " << surf
                                 << " connected to prism; ignoring." << endl;
772 773
                            continue;
                        }
Dave Moxey's avatar
Dave Moxey committed
774

775
                        if (splitEls.count(el->GetId()) > 0)
776
                        {
777 778
                            cerr << "WARNING: prism already found; "
                                 << "ignoring" << endl;
779 780
                        }

781
                        splitEls[el->GetId()] = j;
Dave Moxey's avatar
Dave Moxey committed
782
                    }
783
                    else if (validElTypes.count(el->GetConf().m_e) == 0)
784
                    {
785 786 787
                        cerr << "WARNING: Unsupported element type "
                             << "found in surface " << j << "; "
                             << "ignoring" << endl;
788 789 790
                        continue;
                    }
                }
Dave Moxey's avatar
Dave Moxey committed
791
            }
792 793 794 795 796 797 798 799 800
        }
    }
    else
    {
        // Otherwise, add all prismatic elements and assume face 1 of
        // the prism lies on the surface.
        for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
        {
            ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
801

802 803 804 805 806
            if (el->GetConf().m_e == LibUtilities::ePrism)
            {
                splitEls[el->GetId()] = 1;
            }
            else if (validElTypes.count(el->GetConf().m_e) > 0)
807
            {
808
                splitEls[el->GetId()] = 0;
809
            }
810 811 812 813 814 815
            else
            {
                continue;
            }
        }
    }
816

817 818 819 820 821 822 823 824 825 826
    if (splitEls.size() == 0)
    {
        cerr << "WARNING: No elements detected to split." << endl;
        return;
    }

    // Erase all elements from the element list. Elements will be
    // re-added as they are split.
    vector<ElementSharedPtr> el = m_mesh->m_element[m_mesh->m_expDim];
    m_mesh->m_element[m_mesh->m_expDim].clear();
827

828 829 830 831
    map<int, SpatialDomains::Geometry3DSharedPtr> geomMap;
    for (int i = 0; i < el.size(); ++i)
    {
        const int elId = el[i]->GetId();
832
        if (splitEls.find(elId) == splitEls.end())
833 834 835 836 837
        {
            continue;
        }

        // Get elemental geometry object and put into map.
838
        geomMap[elId] = std::dynamic_pointer_cast<SpatialDomains::Geometry3D>(
839 840 841
            el[i]->GetGeom(m_mesh->m_spaceDim));
    }

842 843 844 845
    // Iterate over list of elements of expansion dimension.
    for (int i = 0; i < el.size(); ++i)
    {
        const int elId = el[i]->GetId();
846
        auto sIt       = splitEls.find(elId);
847 848 849 850 851 852 853

        if (sIt == splitEls.end())
        {
            m_mesh->m_element[m_mesh->m_expDim].push_back(el[i]);
            continue;
        }

854 855
        SpatialDomains::Geometry3DSharedPtr geom = geomMap[elId];

856
        const int faceNum              = sIt->second;
857 858
        LibUtilities::ShapeType elType = el[i]->GetConf().m_e;

859
        SplitMapHelper &sMap   = splitMap[elType][faceNum];
860 861 862 863 864 865 866 867
        SplitEdgeHelper &sEdge = splitEdge[elType][faceNum];

        // Find quadrilateral boundary faces if any
        std::map<int, int> bLink;
        for (int j = 0; j < sMap.bfacesSize; ++j)
        {
            int bl = el[i]->GetBoundaryLink(sMap.bfaces[j]);
            if (bl != -1)
868
            {
869 870 871
                bLink[sMap.bfaces[j]] = bl;
            }
        }
872

873 874
        // Determine whether to use reverse points.
        LibUtilities::PointsType t =
875 876
            revPoints[elType][faceNum] ? LibUtilities::eBoundaryLayerPoints
                                       : LibUtilities::eBoundaryLayerPointsRev;
877

878
        // Determine value of r based on geometry.
879
        if (ratioIsString)
880
        {
881
            NekDouble x, y, z;
882 883
            NekDouble x1, y1, z1;
            int nverts = geom->GetNumVerts();
884

885
            x = y = z = 0.0;
886

887
            for (int i = 0; i < nverts; ++i)
888
            {
889 890 891 892
                geom->GetVertex(i)->GetCoords(x1, y1, z1);
                x += x1;
                y += y1;
                z += z1;
893
            }
894 895 896 897
            x /= (NekDouble)nverts;
            y /= (NekDouble)nverts;
            z /= (NekDouble)nverts;
            r = rEval.Evaluate(rExprId, x, y, z, 0.0);
898
        }
899

900
        LocalRegions::ExpansionSharedPtr q;
901

902 903 904 905
        if (elType == LibUtilities::ePrism)
        {
            // Create basis.
            LibUtilities::BasisKey B0(
906 907 908 909
                LibUtilities::eModified_A, nq, LibUtilities::PointsKey(nq, pt));
            LibUtilities::BasisKey B1(LibUtilities::eModified_A,
                                      2,
                                      LibUtilities::PointsKey(nl + 1, t, r));
910
            LibUtilities::BasisKey B2(
911
                LibUtilities::eModified_B, nq, LibUtilities::PointsKey(nq, pt));
912 913 914

            // Create local region.
            SpatialDomains::PrismGeomSharedPtr g =
915
                std::dynamic_pointer_cast<SpatialDomains::PrismGeom>(geom);
916 917 918 919 920 921 922
            q = MemoryManager<LocalRegions::PrismExp>::AllocateSharedPtr(
                B0, B1, B2, g);
        }
        else if (elType == LibUtilities::eHexahedron)
        {
            // Create basis.
            LibUtilities::BasisKey B0(
923 924 925 926
                LibUtilities::eModified_A, nq, LibUtilities::PointsKey(nq, pt));
            LibUtilities::BasisKey B1(LibUtilities::eModified_A,
                                      2,
                                      LibUtilities::PointsKey(nl + 1, t, r));
927 928 929

            // Create local region.
            SpatialDomains::HexGeomSharedPtr g =
930
                std::dynamic_pointer_cast<SpatialDomains::HexGeom>(geom);
931 932 933
            q = MemoryManager<LocalRegions::HexExp>::AllocateSharedPtr(
                B0, B0, B1, g);
        }
934

935
        // Grab co-ordinates.
936 937 938 939
        Array<OneD, NekDouble> x(nq * nq * (nl + 1));
        Array<OneD, NekDouble> y(nq * nq * (nl + 1));
        Array<OneD, NekDouble> z(nq * nq * (nl + 1));
        q->GetCoords(x, y, z);
940

941 942
        int nSplitEdge = sEdge.size;
        vector<vector<NodeSharedPtr> > edgeNodes(nSplitEdge);
943

944 945 946 947 948
        // Loop over edges to be split.
        for (int j = 0; j < nSplitEdge; ++j)
        {
            int locEdge = sEdge.edge[j];
            int edgeId  = el[i]->GetEdge(locEdge)->m_id;
949

950 951
            // Determine whether we have already generated vertices
            // along this edge.
952
            auto eIt = edgeMap.find(edgeId);
953

954 955 956
            if (eIt == edgeMap.end())
            {
                // If not then resize storage to hold new points.
957
                edgeNodes[j].resize(nl + 1);
958

959 960 961 962
                // Re-use existing vertices at endpoints of edge to
                // avoid duplicating the existing vertices.
                edgeNodes[j][0]  = el[i]->GetVertex(sEdge.edgeVert[j][0]);
                edgeNodes[j][nl] = el[i]->GetVertex(sEdge.edgeVert[j][1]);
963

964
                // Variable geometric ratio
965
                if (ratioIsString)
966
                {
967 968 969
                    NekDouble x0, y0, z0;
                    NekDouble x1, y1, z1;
                    NekDouble xm, ym, zm;
970

971 972 973 974
                    // -> Find edge end and mid points
                    x0 = x[sEdge.offset[j]];
                    y0 = y[sEdge.offset[j]];
                    z0 = z[sEdge.offset[j]];
Dave Moxey's avatar
Dave Moxey committed
975

976 977 978
                    x1 = x[sEdge.offset[j] + nl * nq];
                    y1 = y[sEdge.offset[j] + nl * nq];
                    z1 = z[sEdge.offset[j] + nl * nq];
979

980 981 982
                    xm = 0.5 * (x0 + x1);
                    ym = 0.5 * (y0 + y1);
                    zm = 0.5 * (z0 + z1);
983

984 985
                    // evaluate r factor based on mid point value
                    NekDouble rnew;
986
                    rnew = rEval.Evaluate(rExprId, xm, ym, zm, 0.0);
987

988
                    // Get basis with new r;
989 990 991
                    LibUtilities::PointsKey Pkey(nl + 1, t, rnew);
                    LibUtilities::PointsSharedPtr newP =
                        LibUtilities::PointsManager()[Pkey];
992 993 994 995 996

                    const Array<OneD, const NekDouble> z = newP->GetZ();

                    // Create new interior nodes based on this new blend
                    for (int k = 1; k < nl; ++k)
997
                    {
998 999 1000 1001 1002
                        xm = 0.5 * (1 + z[k]) * (x1 - x0) + x0;
                        ym = 0.5 * (1 + z[k]) * (y1 - y0) + y0;
                        zm = 0.5 * (1 + z[k]) * (z1 - z0) + z0;
                        edgeNodes[j][k] =
                            NodeSharedPtr(new Node(nodeId++, xm, ym, zm));
1003 1004
                    }
                }
1005
                else
1006
                {
1007 1008
                    // Create new interior nodes.
                    for (int k = 1; k < nl; ++k)
1009
                    {
1010
                        int pos         = sEdge.offset[j] + k * sEdge.inc[j];
1011 1012
                        edgeNodes[j][k] = NodeSharedPtr(
                            new Node(nodeId++, x[pos], y[pos], z[pos]));
1013
                    }
1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035
                }

                // Store these edges in edgeMap.
                edgeMap[edgeId] = edgeNodes[j];
            }
            else
            {
                edgeNodes[j] = eIt->second;
            }
        }

        // Create element layers.
        for (int j = 0; j < nl; ++j)
        {
            // Offset of this layer within the collapsed coordinate
            // system.
            int offset = j * sMap.layerOff;

            // Get corner vertices.
            vector<NodeSharedPtr> nodeList(sMap.size);
            for (int k = 0; k < sMap.size; ++k)
            {
1036
                nodeList[k] = edgeNodes[sMap.conn[k][0]][j + sMap.conn[k][1]];
1037 1038 1039 1040
            }

            // Create the element.
            ElmtConfig conf(elType, 1, true, true, false);
1041 1042
            ElementSharedPtr elmt = GetElementFactory().CreateInstance(
                elType, conf, nodeList, el[i]->GetTagList());
1043 1044 1045 1046

            // Add high order nodes to split prismatic edges.
            for (int l = 0; l < sMap.size; ++l)
            {
1047 1048
                EdgeSharedPtr HOedge = elmt->GetEdge(sMap.edge[l]);
                for (int k = 1; k < nq - 1; ++k)
1049
                {
1050 1051 1052
                    int pos = offset + sMap.offset[l] + k * sMap.inc[l];
                    HOedge->m_edgeNodes.push_back(NodeSharedPtr(
                        new Node(nodeId++, x[pos], y[pos], z[pos])));
1053 1054 1055
                }
                HOedge->m_curveType = pt;
            }
1056

1057 1058
            // Change the surface elements to match the layers of
            // elements on the boundary of the domain.
1059
            for (auto &it : bLink)
1060
            {
1061 1062
                int fid = it.first;
                int bl  = it.second;
1063

1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080
                vector<NodeSharedPtr> qNodeList(4);
                for (int k = 0; k < 4; ++k)
                {
                    qNodeList[k] = nodeList[faceNodeMap[elType][fid][k]];
                }
                vector<int> tagBE;
                tagBE =
                    m_mesh->m_element[m_mesh->m_expDim - 1][bl]->GetTagList();
                ElmtConfig bconf(
                    LibUtilities::eQuadrilateral, 1, true, true, false);
                ElementSharedPtr boundaryElmt =
                    GetElementFactory().CreateInstance(
                        LibUtilities::eQuadrilateral, bconf, qNodeList, tagBE);

                // Overwrite first layer boundary element with new
                // boundary element, otherwise push this back to end of
                // the boundary list
1081 1082
                if (j == 0)
                {
1083
                    m_mesh->m_element[m_mesh->m_expDim - 1][bl] = boundaryElmt;
1084 1085 1086
                }
                else
                {
1087 1088
                    m_mesh->m_element[m_mesh->m_expDim - 1].push_back(
                        boundaryElmt);
1089 1090
                }
            }
1091

1092
            m_mesh->m_element[m_mesh->m_expDim].push_back(elmt);
1093 1094
        }
    }
1095 1096
}
}
1097
}