ProcessEquiSpacedOutput.cpp

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////////
//
//  File: ProcessEquiSpacedOutput.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.
//
//  Description: Set up fields as interpolation to equispaced output
//
////////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <string>
using namespace std;

#include "ProcessEquiSpacedOutput.h"

#include <LibUtilities/BasicUtils/ParseUtils.hpp>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Foundations/Interp.h>
#include <StdRegions/StdQuadExp.h>
#include <StdRegions/StdTriExp.h>
#include <boost/math/special_functions/fpclassify.hpp>

namespace Nektar
{
namespace FieldUtils
{

ModuleKey ProcessEquiSpacedOutput::className =
    GetModuleFactory().RegisterCreatorFunction(
        ModuleKey(eProcessModule, "equispacedoutput"),
        ProcessEquiSpacedOutput::create,
        "Write data as equi-spaced output using simplices to represent the "
        "data for connecting points");

ProcessEquiSpacedOutput::ProcessEquiSpacedOutput(FieldSharedPtr f)
    : ProcessModule(f)
{
    f->m_setUpEquiSpacedFields = true;

    m_config["tetonly"] =
        ConfigOption(true, "NotSet", "Only process tetrahedral elements");

    m_config["modalenergy"] =
        ConfigOption(true, "NotSet", "Write output as modal energy");
}

ProcessEquiSpacedOutput::~ProcessEquiSpacedOutput()
{
}

void ProcessEquiSpacedOutput::Process(po::variables_map &vm)
{
    SetupEquiSpacedField();
}

void ProcessEquiSpacedOutput::SetupEquiSpacedField(void)
{
    if (m_f->m_verbose)
    {
        if (m_f->m_comm->TreatAsRankZero())
        {
            cout << "Interpolating fields to equispaced..." << endl;
        }
    }
    int nel = m_f->m_exp[0]->GetExpSize();
    if (!nel)
    {
        m_f->m_fieldPts = LibUtilities::NullPtsField;
        return;
    }

    int coordim  = m_f->m_exp[0]->GetCoordim(0);
    int shapedim = m_f->m_exp[0]->GetExp(0)->GetShapeDimension();
    int npts     = m_f->m_exp[0]->GetTotPoints();
    Array<OneD, Array<OneD, NekDouble> > coords(3);

    // Check if we have a homogeneous expansion
    bool homogeneous1D = false;
    if (m_f->m_fielddef.size())
    {
        if (m_f->m_fielddef[0]->m_numHomogeneousDir == 1)
        {
            coordim++;
            shapedim++;
            homogeneous1D = true;
        }
        else if (m_f->m_fielddef[0]->m_numHomogeneousDir == 2)
        {
            ASSERTL0(false, "Homegeneous2D case not supported");
        }
    }
    else
    {
        if (m_f->m_session->DefinesSolverInfo("HOMOGENEOUS"))
        {
            std::string HomoStr = m_f->m_session->GetSolverInfo("HOMOGENEOUS");

            if ((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D") ||
                (HomoStr == "1D") || (HomoStr == "Homo1D"))
            {
                coordim++;
                shapedim++;
                homogeneous1D = true;
            }
            if ((HomoStr == "HOMOGENEOUS2D") || (HomoStr == "Homogeneous2D") ||
                (HomoStr == "2D") || (HomoStr == "Homo2D"))
            {
                ASSERTL0(false, "Homegeneous2D case not supported");
            }
        }
    }

    // set up the number of points in each element
    int newpoints    = 0;
    int newtotpoints = 0;

    Array<OneD, int> conn;
    int prevNcoeffs = 0;
    int prevNpoints = 0;
    int cnt         = 0;

    // identify face 1 connectivity for prisms
    map<int, StdRegions::Orientation> face0orient;
    set<int> prismorient;
    LocalRegions::ExpansionSharedPtr e;

    // prepare PtsField
    vector<std::string> fieldNames;
    vector<int> ppe;
    vector<Array<OneD, int> > ptsConn;
    int nfields;

    for (int i = 0; i < nel; ++i)
    {
        e = m_f->m_exp[0]->GetExp(i);
        if (e->DetShapeType() == LibUtilities::ePrism)
        {
            StdRegions::Orientation forient = e->GetForient(0);
            int fid                         = e->GetGeom()->GetFid(0);
            if (face0orient.count(fid))
            { // face 1 meeting face 1 so reverse this id
                prismorient.insert(i);
            }
            else
            {
                // just store if Dir 1 is fwd or bwd
                if ((forient == StdRegions::eDir1BwdDir1_Dir2FwdDir2) ||
                    (forient == StdRegions::eDir1BwdDir1_Dir2BwdDir2) ||
                    (forient == StdRegions::eDir1BwdDir2_Dir2FwdDir1) ||
                    (forient == StdRegions::eDir1BwdDir2_Dir2BwdDir1))
                {
                    face0orient[fid] = StdRegions::eBwd;
                }
                else
                {
                    face0orient[fid] = StdRegions::eFwd;
                }
            }
        }
    }

    for (int i = 0; i < nel; ++i)
    {
        e = m_f->m_exp[0]->GetExp(i);
        if (e->DetShapeType() == LibUtilities::ePrism)
        {
            int fid = e->GetGeom()->GetFid(2);
            // check to see if face 2 meets face 1
            if (face0orient.count(fid))
            {
                // check to see how face 2 is orientated
                StdRegions::Orientation forient2 = e->GetForient(2);
                StdRegions::Orientation forient0 = face0orient[fid];

                // If dir 1 or forient2 is bwd then check agains
                // face 1 value
                if ((forient2 == StdRegions::eDir1BwdDir1_Dir2FwdDir2) ||
                    (forient2 == StdRegions::eDir1BwdDir1_Dir2BwdDir2) ||
                    (forient2 == StdRegions::eDir1BwdDir2_Dir2FwdDir1) ||
                    (forient2 == StdRegions::eDir1BwdDir2_Dir2BwdDir1))
                {
                    if (forient0 == StdRegions::eFwd)
                    {
                        prismorient.insert(i);
                    }
                }
                else
                {
                    if (forient0 == StdRegions::eBwd)
                    {
                        prismorient.insert(i);
                    }
                }
            }
        }
    }

    for (int i = 0; i < nel; ++i)
    {
        e = m_f->m_exp[0]->GetExp(i);
        if (m_config["tetonly"].m_beenSet)
        {
            if (m_f->m_exp[0]->GetExp(i)->DetShapeType() !=
                LibUtilities::eTetrahedron)
            {
                continue;
            }
        }

        switch (e->DetShapeType())
        {
            case LibUtilities::eSegment:
            {
                int npoints0 = e->GetBasis(0)->GetNumPoints();

                newpoints =
                    LibUtilities::StdSegData::getNumberOfCoefficients(npoints0);
            }
            break;
            case LibUtilities::eTriangle:
            {
                int np0 = e->GetBasis(0)->GetNumPoints();
                int np1 = e->GetBasis(1)->GetNumPoints();
                int np  = max(np0, np1);
                newpoints =
                    LibUtilities::StdTriData::getNumberOfCoefficients(np, np);
            }
            break;
            case LibUtilities::eQuadrilateral:
            {
                int np0 = e->GetBasis(0)->GetNumPoints();
                int np1 = e->GetBasis(1)->GetNumPoints();
                int np  = max(np0, np1);

                newpoints =
                    LibUtilities::StdQuadData::getNumberOfCoefficients(np, np);
            }
            break;
            case LibUtilities::eTetrahedron:
            {
                int np0 = e->GetBasis(0)->GetNumPoints();
                int np1 = e->GetBasis(1)->GetNumPoints();
                int np2 = e->GetBasis(2)->GetNumPoints();
                int np  = max(np0, max(np1, np2));

                newpoints = LibUtilities::StdTetData::getNumberOfCoefficients(
                    np, np, np);
            }
            break;
            case LibUtilities::ePrism:
            {
                int np0 = e->GetBasis(0)->GetNumPoints();
                int np1 = e->GetBasis(1)->GetNumPoints();
                int np2 = e->GetBasis(2)->GetNumPoints();
                int np  = max(np0, max(np1, np2));

                newpoints = LibUtilities::StdPrismData::getNumberOfCoefficients(
                    np, np, np);
            }
            break;
            case LibUtilities::ePyramid:
            {
                int np0 = e->GetBasis(0)->GetNumPoints();
                int np1 = e->GetBasis(1)->GetNumPoints();
                int np2 = e->GetBasis(2)->GetNumPoints();
                int np  = max(np0, max(np1, np2));

                newpoints = LibUtilities::StdPyrData::getNumberOfCoefficients(
                    np, np, np);
            }
            break;
            case LibUtilities::eHexahedron:
            {
                int np0 = e->GetBasis(0)->GetNumPoints();
                int np1 = e->GetBasis(1)->GetNumPoints();
                int np2 = e->GetBasis(2)->GetNumPoints();
                int np  = max(np0, max(np1, np2));

                newpoints = LibUtilities::StdHexData::getNumberOfCoefficients(
                    np, np, np);
            }
            break;
            default:
            {
                ASSERTL0(false, "Points not known");
            }
        }

        ppe.push_back(newpoints);
        newtotpoints += newpoints;

        if (e->DetShapeType() == LibUtilities::ePrism)
        {
            bool standard = true;

            if (prismorient.count(i))
            {
                standard = false; // reverse direction
            }

            e->GetSimplexEquiSpacedConnectivity(conn, standard);
        }
        else
        {

            if ((prevNcoeffs != e->GetNcoeffs()) ||
                (prevNpoints != e->GetTotPoints()))
            {
                prevNcoeffs = e->GetNcoeffs();
                prevNpoints = e->GetTotPoints();

                e->GetSimplexEquiSpacedConnectivity(conn);
            }
        }
        Array<OneD, int> newconn(conn.num_elements());
        for (int j = 0; j < conn.num_elements(); ++j)
        {
            newconn[j] = conn[j] + cnt;
        }

        ptsConn.push_back(newconn);
        cnt += newpoints;
    }

    if (m_f->m_fielddef.size())
    {
        nfields = m_f->m_exp.size();
    }
    else // just the mesh points
    {
        nfields = 0;
    }

    Array<OneD, Array<OneD, NekDouble> > pts(nfields + coordim);

    for (int i = 0; i < nfields + coordim; ++i)
    {
        pts[i] = Array<OneD, NekDouble>(newtotpoints);
    }

    // Interpolate coordinates
    for (int i = 0; i < coordim; ++i)
    {
        coords[i] = Array<OneD, NekDouble>(npts);
    }

    for (int i = coordim; i < 3; ++i)
    {
        coords[i] = NullNekDouble1DArray;
    }

    m_f->m_exp[0]->GetCoords(coords[0], coords[1], coords[2]);

    int nq1 = m_f->m_exp[0]->GetTotPoints();

    Array<OneD, NekDouble> x1(nq1);
    Array<OneD, NekDouble> y1(nq1);
    Array<OneD, NekDouble> z1(nq1);

    m_f->m_exp[0]->GetCoords(x1, y1, z1);

    Array<OneD, NekDouble> tmp;

    for (int n = 0; n < coordim; ++n)
    {
        cnt      = 0;
        int cnt1 = 0;
        for (int i = 0; i < nel; ++i)
        {
            m_f->m_exp[0]->GetExp(i)->PhysInterpToSimplexEquiSpaced(
                coords[n] + cnt, tmp = pts[n] + cnt1);
            cnt1 += ppe[i];
            cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
        }
    }

    if (m_f->m_fielddef.size())
    {
        ASSERTL0(m_f->m_fielddef[0]->m_fields.size() == m_f->m_exp.size(),
                 "More expansion defined than fields");

        for (int n = 0; n < m_f->m_exp.size(); ++n)
        {
            cnt      = 0;
            int cnt1 = 0;

            if (m_config["modalenergy"].m_beenSet)
            {
                Array<OneD, const NekDouble> phys = m_f->m_exp[n]->GetPhys();
                for (int i = 0; i < nel; ++i)
                {
                    GenOrthoModes(i, phys + cnt, tmp = pts[coordim + n] + cnt1);
                    cnt1 += ppe[i];
                    cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
                }
            }
            else
            {
                Array<OneD, const NekDouble> phys = m_f->m_exp[n]->GetPhys();
                for (int i = 0; i < nel; ++i)
                {
                    m_f->m_exp[0]->GetExp(i)->PhysInterpToSimplexEquiSpaced(
                        phys + cnt, tmp = pts[coordim + n] + cnt1);
                    cnt1 += ppe[i];
                    cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
                }
            }

            // Set up Variable string.
            fieldNames.push_back(m_f->m_fielddef[0]->m_fields[n]);
        }
    }

    m_f->m_fieldPts = MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(
        coordim, fieldNames, pts);
    if (shapedim == 1)
    {
        m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsSegBlock);
    }

    if (shapedim == 2)
    {
        m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTriBlock);
    }
    else if (shapedim == 3)
    {
        m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTetBlock);
    }
    m_f->m_fieldPts->SetConnectivity(ptsConn);
    if (homogeneous1D)
    {
        SetHomogeneousConnectivity();
    }
}

void ProcessEquiSpacedOutput::SetHomogeneousConnectivity(void)
{
    LocalRegions::ExpansionSharedPtr e;
    int nel          = m_f->m_exp[0]->GetPlane(0)->GetExpSize();
    int nPlanes      = m_f->m_exp[0]->GetZIDs().num_elements();
    int npts         = m_f->m_fieldPts->GetNpoints();
    int nptsPerPlane = npts / nPlanes;
    int coordim      = 3;

    if (m_f->m_exp[0]->GetHomogeneousBasis()->GetBasisType() ==
            LibUtilities::eFourier &&
        m_f->m_exp[0]->GetHomogeneousBasis()->GetPointsType() ==
            LibUtilities::eFourierEvenlySpaced)
    {
        // Write points with extra plane
        Array<OneD, Array<OneD, NekDouble> > pts;
        m_f->m_fieldPts->GetPts(pts);
        Array<OneD, Array<OneD, NekDouble> > newPts(pts.num_elements());
        for (int i = 0; i < pts.num_elements(); i++)
        {
            newPts[i] = Array<OneD, NekDouble>(npts + nptsPerPlane);
            // Copy old points
            Vmath::Vcopy(npts, pts[i], 1, newPts[i], 1);
            // Get new plane
            Array<OneD, NekDouble> extraPlane(nptsPerPlane, newPts[i] + npts);
            if (m_f->m_comm->GetColumnComm()->GetSize() == 1)
            {
                if ( i == (coordim-1))
                {
                    // z-coordinate
                    Vmath::Sadd(nptsPerPlane, m_f->m_exp[0]->GetHomoLen(),
                                pts[i], 1, extraPlane, 1);
                }
                else
                {
                    Vmath::Vcopy(nptsPerPlane, pts[i], 1, extraPlane, 1);
                }
            }
            else
            {
                // Determine to and from rank for communication
                int size     = m_f->m_comm->GetColumnComm()->GetSize();
                int rank     = m_f->m_comm->GetColumnComm()->GetRank();
                int fromRank = (rank + 1) % size;
                int toRank   = (rank == 0) ? size - 1 : rank - 1;
                // Communicate using SendRecv
                Array<OneD, NekDouble> send(nptsPerPlane, pts[i]);
                m_f->m_comm->GetColumnComm()->SendRecv(toRank, send, fromRank,
                                                       extraPlane);
                // Adjust z-coordinate of last process
                if (i == (coordim-1) && (rank == size - 1) )
                {
                    Vmath::Sadd(nptsPerPlane, m_f->m_exp[0]->GetHomoLen(),
                                extraPlane, 1, extraPlane, 1);
                }
            }
        }
        m_f->m_fieldPts->SetPts(newPts);
        // Now extend plane connectivity
        vector<Array<OneD, int> > oldConn;
        Array<OneD, int> conn;
        m_f->m_fieldPts->GetConnectivity(oldConn);
        vector<Array<OneD, int> > newConn = oldConn;
        int connPerPlane = oldConn.size() / nPlanes;
        for (int i = 0; i < connPerPlane; i++)
        {
            conn = Array<OneD, int>(oldConn[i].num_elements());
            for (int j = 0; j < conn.num_elements(); j++)
            {
                conn[j] = oldConn[i][j] + npts;
            }
            newConn.push_back(conn);
        }
        m_f->m_fieldPts->SetConnectivity(newConn);

        nPlanes++;
        npts += nptsPerPlane;
    }

    vector<Array<OneD, int> > oldConn;
    vector<Array<OneD, int> > newConn;
    Array<OneD, int>          conn;
    m_f->m_fieldPts->GetConnectivity(oldConn);

    // 2DH1D case
    if (m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTriBlock)
    {
        for(int i = 0; i < nel; ++i)
        {
            int nLines = oldConn[i].num_elements()/2;
            // Create array for new connectivity
            // (2 triangles between each plane for each line)
            conn = Array<OneD, int> (2*3*nLines*(nPlanes-1));
            int cnt = 0;
            for(int n = 0; n<nLines; n++)
            {
                // Define new connectivity with triangles
                for ( int p = 0; p<nPlanes-1; p++)
                {
                    conn[cnt++] = oldConn[i+  p  *nel][2*n+0];
                    conn[cnt++] = oldConn[i+  p  *nel][2*n+1];
                    conn[cnt++] = oldConn[i+(p+1)*nel][2*n+0];

                    conn[cnt++] = oldConn[i+  p  *nel][2*n+1];
                    conn[cnt++] = oldConn[i+(p+1)*nel][2*n+0];
                    conn[cnt++] = oldConn[i+(p+1)*nel][2*n+1];
                }
            }
            newConn.push_back(conn);
        }
    }
    else if(m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTetBlock)
    {
        // Get maximum number of points per direction in the mesh
        int maxN = 0;
        for(int i = 0; i < nel; ++i)
        {
            e = m_f->m_exp[0]->GetPlane(0)->GetExp(i);

            int np0 = e->GetBasis(0)->GetNumPoints();
            int np1 = e->GetBasis(1)->GetNumPoints();
            int np = max(np0,np1);
            maxN = max(np, maxN);
        }

        // Create a global numbering for points in plane 0
        Array<OneD, int> vId(nptsPerPlane);
        int cnt1=0;
        int cnt2=0;
        for(int i = 0; i < nel; ++i)
        {
            e = m_f->m_exp[0]->GetPlane(0)->GetExp(i);
            int np0 = e->GetBasis(0)->GetNumPoints();
            int np1 = e->GetBasis(1)->GetNumPoints();
            int np = max(np0,np1);
            switch(e->DetShapeType())
            {
            case LibUtilities::eTriangle:
                {
                    int newpoints = LibUtilities::StdTriData::
                                        getNumberOfCoefficients(np,np);

                    // Vertex numbering
                    vId[cnt1]                 = e->GetGeom()->GetVid(0);
                    vId[cnt1+np-1]            = e->GetGeom()->GetVid(1);
                    vId[cnt1+newpoints-1]     = e->GetGeom()->GetVid(2);

                    // Edge numbering
                    StdRegions::Orientation             edgeOrient;
                    int edge1 = 0;
                    int edge2 = 0;
                    for (int n = 1; n < np-1; n++)
                    {
                        // Edge 0
                        edgeOrient          = e->GetGeom()->GetEorient(0);
                        if (edgeOrient==StdRegions::eForwards)
                        {
                            vId[cnt1+n] = 4*nel +
                                            maxN*e->GetGeom()->GetEid(0) + n;
                        }
                        else
                        {
                            vId[cnt1+np-1-n] = 4*nel +
                                            maxN*e->GetGeom()->GetEid(0) + n;
                        }

                        // Edge 1
                        edgeOrient          = e->GetGeom()->GetEorient(1);
                        if (edgeOrient==StdRegions::eForwards)
                        {
                            edge1 += np-n;
                            vId[cnt1+np-1+edge1] = 4*nel +
                                            maxN*e->GetGeom()->GetEid(1) + n;
                        }
                        else
                        {
                            edge1 += n;
                            vId[cnt1+newpoints-1-edge1] = 4*nel +
                                            maxN*e->GetGeom()->GetEid(1) + n;
                        }

                        // Edge 2
                        edgeOrient          = e->GetGeom()->GetEorient(2);
                        if (edgeOrient==StdRegions::eForwards)
                        {
                            edge2 += n+1;
                            vId[cnt1+newpoints-1-edge2] = 4*nel +
                                            maxN*e->GetGeom()->GetEid(2) + n;
                        }
                        else
                        {
                            edge2 += np+1-n;
                            vId[cnt1+edge2] = 4*nel +
                                            maxN*e->GetGeom()->GetEid(2) + n;
                        }
                    }

                    // Interior numbering
                    edge2 = 0;
                    for (int n = 1; n < np-1; n++)
                    {
                        edge2 += np+1-n;
                        for (int m = 1; m < np-n-1; m++)
                        {
                            vId[cnt1+edge2+m] = 4*nel + maxN*4*nel + cnt2;
                            cnt2++;
                        }
                    }
                    cnt1+= newpoints;
                }
                break;
            case LibUtilities::eQuadrilateral:
                {
                    int newpoints = LibUtilities::StdQuadData::
                                        getNumberOfCoefficients(np,np);
                    // Vertex numbering
                    vId[cnt1]           = e->GetGeom()->GetVid(0);
                    vId[cnt1+np-1]      = e->GetGeom()->GetVid(1);
                    vId[cnt1+np*np-1]   = e->GetGeom()->GetVid(2);
                    vId[cnt1+np*(np-1)] = e->GetGeom()->GetVid(3);

                    // Edge numbering
                    StdRegions::Orientation             edgeOrient;
                    for (int n = 1; n < np-1; n++)
                    {
                        // Edge 0
                        edgeOrient          = e->GetGeom()->GetEorient(0);
                        if (edgeOrient==StdRegions::eForwards)
                        {
                            vId[cnt1+n] = 4*nel +
                                            maxN*e->GetGeom()->GetEid(0) + n;
                        }
                        else
                        {
                            vId[cnt1+np-1-n] = 4*nel +
                                            maxN*e->GetGeom()->GetEid(0) + n;
                        }

                        // Edge 1
                        edgeOrient          = e->GetGeom()->GetEorient(1);
                        if (edgeOrient==StdRegions::eForwards)
                        {
                            vId[cnt1+np-1+n*np] = 4*nel +
                                            maxN*e->GetGeom()->GetEid(1) + n;
                        }
                        else
                        {
                            vId[cnt1+np*np-1-n*np] = 4*nel +
                                            maxN*e->GetGeom()->GetEid(1) + n;
                        }

                        // Edge 2
                        edgeOrient          = e->GetGeom()->GetEorient(2);
                        if (edgeOrient==StdRegions::eForwards)
                        {
                            vId[cnt1+np*np-1-n] = 4*nel +
                                            maxN*e->GetGeom()->GetEid(2) + n;
                        }
                        else
                        {
                            vId[cnt1+np*(np-1)+n] = 4*nel +
                                            maxN*e->GetGeom()->GetEid(2) + n;
                        }

                        // Edge 3
                        edgeOrient          = e->GetGeom()->GetEorient(3);
                        if (edgeOrient==StdRegions::eForwards)
                        {
                            vId[cnt1+np*(np-1)-n*np] = 4*nel +
                                            maxN*e->GetGeom()->GetEid(3) + n;
                        }
                        else
                        {
                            vId[cnt1+n*np] = 4*nel +
                                            maxN*e->GetGeom()->GetEid(3) + n;
                        }
                    }

                    // Interior numbering
                    for (int n = 1; n < np-1; n++)
                    {
                        for (int m = 1; m < np-1; m++)
                        {
                            vId[cnt1+m*np+n] = 4*nel + maxN*4*nel + cnt2;
                            cnt2++;
                        }
                    }
                    cnt1+= newpoints;
                }
                break;
            default:
                {
                    ASSERTL0(false,"Points not known");
                }
            }
        }

        // Crete new connectivity using homogeneous information
        for(int i = 0; i < nel; ++i)
        {
            int nTris = oldConn[i].num_elements()/3;
            // Create array for new connectivity
            // (3 tetrahedra between each plane for each triangle)
            conn = Array<OneD, int> (4*3*nTris*(nPlanes-1));
            cnt2=0;
            for(int n = 0; n<nTris; n++)
            {
                // Get id of vertices in this triangle (on plane 0)
                int vId0 = vId[oldConn[i][3*n+0]];
                int vId1 = vId[oldConn[i][3*n+1]];
                int vId2 = vId[oldConn[i][3*n+2]];

                // Determine ordering based on global Id of pts
                Array<OneD, int> rot(3);
                if ( (vId0<vId1) && (vId0<vId2))
                {
                    rot[0] = 0;
                    if (vId1<vId2)
                    {
                        rot[1] = 1;
                        rot[2] = 2;
                    }
                    else
                    {
                        rot[1] = 2;
                        rot[2] = 1;
                    }
                }
                else if ( (vId1<vId0) && (vId1<vId2))
                {
                    rot[0] = 1;
                    if (vId0<vId2)
                    {
                        rot[1] = 0;
                        rot[2] = 2;
                    }
                    else
                    {
                        rot[1] = 2;
                        rot[2] = 0;
                    }
                }
                else
                {
                    rot[0] = 2;
                    if (vId0<vId1)
                    {
                        rot[1] = 0;
                        rot[2] = 1;
                    }
                    else
                    {
                        rot[1] = 1;
                        rot[2] = 0;
                    }
                }

                // Define new connectivity with tetrahedra
                for ( int p = 0; p<nPlanes-1; p++)
                {
                    conn[cnt2++] = oldConn[i+  p  *nel][3*n+rot[0]];
                    conn[cnt2++] = oldConn[i+  p  *nel][3*n+rot[1]];
                    conn[cnt2++] = oldConn[i+  p  *nel][3*n+rot[2]];
                    conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];

                    conn[cnt2++] = oldConn[i+  p  *nel][3*n+rot[0]];
                    conn[cnt2++] = oldConn[i+  p  *nel][3*n+rot[1]];
                    conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
                    conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[1]];

                    conn[cnt2++] = oldConn[i+  p  *nel][3*n+rot[0]];
                    conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[1]];
                    conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[2]];
                    conn[cnt2++] = oldConn[i+(p+1)*nel][3*n+rot[0]];
                }
            }
            newConn.push_back(conn);
        }
    }
    else
    {
        ASSERTL0( false, "Points type not supported.");
    }

    m_f->m_fieldPts->SetConnectivity(newConn);
}

void ProcessEquiSpacedOutput::GenOrthoModes(
    int n,
    const Array<OneD, const NekDouble> &phys,
    Array<OneD, NekDouble> &coeffs)
{
    LocalRegions::ExpansionSharedPtr e;
    e = m_f->m_exp[0]->GetExp(n);

    switch (e->DetShapeType())
    {
        case LibUtilities::eTriangle:
        {
            int np0 = e->GetBasis(0)->GetNumPoints();
            int np1 = e->GetBasis(1)->GetNumPoints();
            int np  = max(np0, np1);

            // to ensure points are correctly projected to np need
            // to increase the order slightly of coordinates
            LibUtilities::PointsKey pa(np + 1, e->GetPointsType(0));
            LibUtilities::PointsKey pb(np, e->GetPointsType(1));
            Array<OneD, NekDouble> tophys(np * (np + 1));

            LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, np, pa);
            LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B, np, pb);
            StdRegions::StdTriExp OrthoExp(Ba, Bb);

            // interpolate points to new phys points!
            LibUtilities::Interp2D(e->GetBasis(0)->GetBasisKey(),
                                   e->GetBasis(1)->GetBasisKey(), phys, Ba, Bb,
                                   tophys);

            OrthoExp.FwdTrans(tophys, coeffs);
            break;
        }
        case LibUtilities::eQuadrilateral:
        {
            int np0 = e->GetBasis(0)->GetNumPoints();
            int np1 = e->GetBasis(1)->GetNumPoints();
            int np  = max(np0, np1);

            LibUtilities::PointsKey pa(np + 1, e->GetPointsType(0));
            LibUtilities::PointsKey pb(np + 1, e->GetPointsType(1));
            Array<OneD, NekDouble> tophys((np + 1) * (np + 1));

            LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, np, pa);
            LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, np, pb);
            StdRegions::StdQuadExp OrthoExp(Ba, Bb);

            // interpolate points to new phys points!
            LibUtilities::Interp2D(e->GetBasis(0)->GetBasisKey(),
                                   e->GetBasis(1)->GetBasisKey(), phys, Ba, Bb,
                                   tophys);

            OrthoExp.FwdTrans(phys, coeffs);
            break;
        }
        default:
            ASSERTL0(false, "Shape needs setting up");
            break;
    }
}
}
}