Nektar::FieldUtils::ProcessEquiSpacedOutput Class Reference

This processing module interpolates one field to another. More...

#include <ProcessEquiSpacedOutput.h>

Inheritance diagram for Nektar::FieldUtils::ProcessEquiSpacedOutput:

Public Member Functions
	ProcessEquiSpacedOutput (FieldSharedPtr f)
	~ProcessEquiSpacedOutput () override
Public Member Functions inherited from Nektar::FieldUtils::ProcessModule
	ProcessModule ()
	ProcessModule (FieldSharedPtr p_f)
Public Member Functions inherited from Nektar::FieldUtils::Module
FIELD_UTILS_EXPORT	Module (FieldSharedPtr p_f)
virtual	~Module ()=default
void	Process (po::variables_map &vm)
std::string	GetModuleName ()
std::string	GetModuleDescription ()
const ConfigOption &	GetConfigOption (const std::string &key) const
ModulePriority	GetModulePriority ()
std::vector< ModuleKey >	GetModulePrerequisites ()
FIELD_UTILS_EXPORT void	RegisterConfig (std::string key, std::string value="")
	Register a configuration option with a module. More...
FIELD_UTILS_EXPORT void	PrintConfig ()
	Print out all configuration options for a module. More...
FIELD_UTILS_EXPORT void	SetDefaults ()
	Sets default configuration options for those which have not been set. More...
FIELD_UTILS_EXPORT void	AddFile (std::string fileType, std::string fileName)
FIELD_UTILS_EXPORT void	EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)

Static Public Member Functions
static std::shared_ptr< Module >	create (FieldSharedPtr f)
	Creates an instance of this class. More...

Static Public Attributes
static ModuleKey	className

Protected Member Functions
void	v_Process (po::variables_map &vm) override
	Write mesh to output file. More...
std::string	v_GetModuleName () override
std::string	v_GetModuleDescription () override
ModulePriority	v_GetModulePriority () override
void	SetHomogeneousConnectivity (void)
void	GenOrthoModes (int n, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &coeffs)
Protected Member Functions inherited from Nektar::FieldUtils::Module
	Module ()
virtual void	v_Process (po::variables_map &vm)
virtual std::string	v_GetModuleName ()
virtual std::string	v_GetModuleDescription ()
virtual ModulePriority	v_GetModulePriority ()
virtual std::vector< ModuleKey >	v_GetModulePrerequisites ()

Additional Inherited Members
Public Attributes inherited from Nektar::FieldUtils::Module
FieldSharedPtr	m_f
	Field object. More...
Protected Attributes inherited from Nektar::FieldUtils::Module
std::map< std::string, ConfigOption >	m_config
	List of configuration values. More...
std::set< std::string >	m_allowedFiles
	List of allowed file formats. More...

Detailed Description

This processing module interpolates one field to another.

Definition at line 46 of file ProcessEquiSpacedOutput.h.

Constructor & Destructor Documentation

ProcessEquiSpacedOutput()

Nektar::FieldUtils::ProcessEquiSpacedOutput::ProcessEquiSpacedOutput

(

FieldSharedPtr

f

)

Definition at line 55 of file ProcessEquiSpacedOutput.cpp.

    : ProcessModule(f)
{
    m_config["tetonly"] =
        ConfigOption(true, "0", "Only process tetrahedral elements");
 
    m_config["modalenergy"] =
        ConfigOption(true, "0", "Write output as modal energy");
}

References Nektar::FieldUtils::Module::m_config.

~ProcessEquiSpacedOutput()

Nektar::FieldUtils::ProcessEquiSpacedOutput::~ProcessEquiSpacedOutput ( )

override

Definition at line 65 of file ProcessEquiSpacedOutput.cpp.

{
}

Member Function Documentation

create()

static std::shared_ptr< Module > Nektar::FieldUtils::ProcessEquiSpacedOutput::create ( FieldSharedPtr  f )

inlinestatic

Creates an instance of this class.

Definition at line 50 of file ProcessEquiSpacedOutput.h.

    {
        return MemoryManager<ProcessEquiSpacedOutput>::AllocateSharedPtr(f);
    }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

GenOrthoModes()

void Nektar::FieldUtils::ProcessEquiSpacedOutput::GenOrthoModes ( int  n, const Array< OneD, const NekDouble > &  phys, Array< OneD, NekDouble > &  coeffs  )

protected

Definition at line 809 of file ProcessEquiSpacedOutput.cpp.

{
    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;
    }
}

References ASSERTL0, Nektar::LibUtilities::eOrtho_A, Nektar::LibUtilities::eOrtho_B, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::StdRegions::StdExpansion::FwdTrans(), Nektar::LibUtilities::Interp2D(), and Nektar::FieldUtils::Module::m_f.

Referenced by v_Process().

SetHomogeneousConnectivity()

void Nektar::FieldUtils::ProcessEquiSpacedOutput::SetHomogeneousConnectivity ( void  )

protected

Definition at line 418 of file ProcessEquiSpacedOutput.cpp.

{
    LocalRegions::ExpansionSharedPtr e;
    int nel          = m_f->m_exp[0]->GetPlane(0)->GetExpSize();
    int nPlanes      = m_f->m_exp[0]->GetZIDs().size();
    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.size());
        for (int i = 0; i < pts.size(); 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].size());
            for (int j = 0; j < conn.size(); 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].size() / 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].size() / 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);
}

References ASSERTL0, Nektar::StdRegions::eForwards, Nektar::LibUtilities::eFourier, Nektar::LibUtilities::eFourierEvenlySpaced, Nektar::LibUtilities::ePtsTetBlock, Nektar::LibUtilities::ePtsTriBlock, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::LibUtilities::StdTriData::getNumberOfCoefficients(), Nektar::LibUtilities::StdQuadData::getNumberOfCoefficients(), Nektar::FieldUtils::Module::m_f, CellMLToNektar.cellml_metadata::p, Vmath::Sadd(), and Vmath::Vcopy().

Referenced by v_Process().

v_GetModuleDescription()

std::string Nektar::FieldUtils::ProcessEquiSpacedOutput::v_GetModuleDescription ( )

inlineoverrideprotectedvirtual

Reimplemented from Nektar::FieldUtils::Module.

Definition at line 68 of file ProcessEquiSpacedOutput.h.

    {
        return "Interpolating fields to equispaced";
    }

v_GetModuleName()

std::string Nektar::FieldUtils::ProcessEquiSpacedOutput::v_GetModuleName ( )

inlineoverrideprotectedvirtual

Reimplemented from Nektar::FieldUtils::Module.

Definition at line 63 of file ProcessEquiSpacedOutput.h.

    {
        return "ProcessEquiSpacedOutput";
    }

v_GetModulePriority()

ModulePriority Nektar::FieldUtils::ProcessEquiSpacedOutput::v_GetModulePriority ( )

inlineoverrideprotectedvirtual

Reimplemented from Nektar::FieldUtils::Module.

Definition at line 73 of file ProcessEquiSpacedOutput.h.

    {
        return eConvertExpToPts;
    }

References Nektar::FieldUtils::eConvertExpToPts.

v_Process()

void Nektar::FieldUtils::ProcessEquiSpacedOutput::v_Process ( po::variables_map &  vm )

overrideprotectedvirtual

Write mesh to output file.

Reimplemented from Nektar::FieldUtils::Module.

Definition at line 69 of file ProcessEquiSpacedOutput.cpp.

{
    m_f->SetUpExp(vm);
 
    int nel = m_f->m_exp[0]->GetExpSize();
    if (!nel)
    {
        // Create empty PtsField
        int nfields = m_f->m_variables.size();
        int coordim = 3;
 
        Array<OneD, Array<OneD, NekDouble>> pts(nfields + coordim);
        for (int i = 0; i < nfields + coordim; ++i)
        {
            pts[i] = Array<OneD, NekDouble>(0);
        }
        vector<Array<OneD, int>> ptsConn;
 
        m_f->m_fieldPts =
            MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(
                coordim, m_f->m_variables, pts);
        m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTetBlock);
        m_f->m_fieldPts->SetConnectivity(ptsConn);
        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_numHomogeneousDir == 1)
    {
        coordim++;
        shapedim++;
        homogeneous1D = true;
    }
    else if (m_f->m_numHomogeneousDir == 2)
    {
        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<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->GetTraceOrient(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->GetTraceOrient(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 && m_config["tetonly"].as<bool>())
        {
            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.size());
        for (int j = 0; j < conn.size(); ++j)
        {
            newconn[j] = conn[j] + cnt;
        }
 
        ptsConn.push_back(newconn);
        cnt += newpoints;
    }
 
    nfields = m_f->m_variables.size();
 
    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]);
 
    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();
        }
    }
 
    for (int n = 0; n < m_f->m_variables.size(); ++n)
    {
        cnt      = 0;
        int cnt1 = 0;
 
        if (m_config["modalenergy"].m_beenSet &&
            m_config["modalenergy"].as<bool>())
        {
            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();
            }
        }
    }
 
    m_f->m_fieldPts = MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(
        coordim, m_f->m_variables, 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();
    }
 
    // Save points per element in the point field
    m_f->m_fieldPts->SetPointsPerElement(ppe);
 
    // Clear m_exp
    m_f->m_exp = vector<MultiRegions::ExpListSharedPtr>();
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::StdRegions::eBwd, Nektar::StdRegions::eDir1BwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1BwdDir1_Dir2FwdDir2, Nektar::StdRegions::eDir1BwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1BwdDir2_Dir2FwdDir1, Nektar::StdRegions::eFwd, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::ePtsSegBlock, Nektar::LibUtilities::ePtsTetBlock, Nektar::LibUtilities::ePtsTriBlock, Nektar::LibUtilities::ePyramid, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eSegment, Nektar::LibUtilities::eTetrahedron, Nektar::LibUtilities::eTriangle, GenOrthoModes(), Nektar::LibUtilities::StdSegData::getNumberOfCoefficients(), Nektar::LibUtilities::StdTriData::getNumberOfCoefficients(), Nektar::LibUtilities::StdQuadData::getNumberOfCoefficients(), Nektar::LibUtilities::StdHexData::getNumberOfCoefficients(), Nektar::LibUtilities::StdTetData::getNumberOfCoefficients(), Nektar::LibUtilities::StdPyrData::getNumberOfCoefficients(), Nektar::LibUtilities::StdPrismData::getNumberOfCoefficients(), Nektar::FieldUtils::Module::m_config, Nektar::FieldUtils::Module::m_f, Nektar::NullNekDouble1DArray, and SetHomogeneousConnectivity().

Member Data Documentation

className

ModuleKey Nektar::FieldUtils::ProcessEquiSpacedOutput::className

static

Initial value:

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

Definition at line 54 of file ProcessEquiSpacedOutput.h.