Converter from fld to vtk. More...

#include <OutputVtk.h>

Inheritance diagram for Nektar::FieldUtils::OutputVtk:

Public Member Functions
	OutputVtk (FieldSharedPtr f)

	~OutputVtk () final=default

Public Member Functions inherited from Nektar::FieldUtils::OutputVtkBase
	OutputVtkBase (FieldSharedPtr f)

virtual	~OutputVtkBase ()

Public Member Functions inherited from Nektar::FieldUtils::OutputFileBase
	OutputFileBase (FieldSharedPtr f)

virtual	~OutputFileBase ()

Public Member Functions inherited from Nektar::FieldUtils::OutputModule
	OutputModule (FieldSharedPtr p_f)

FIELD_UTILS_EXPORT void	OpenStream ()
	Open a file for output. More...

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 ()

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 (const FieldSharedPtr &f)
	Creates an instance of this class. More...

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

Static Public Attributes
static ModuleKey	m_className

Static Public Attributes inherited from Nektar::FieldUtils::OutputVtkBase
static ModuleKey	m_className

Protected Member Functions
virtual std::string	v_GetModuleName () override final

virtual void	v_OutputFromPts (po::variables_map &vm) override final
	Write from pts to output file. More...

virtual void	v_OutputFromExp (po::variables_map &vm) override final
	Write from m_exp to output file. More...

virtual void	v_OutputFromData (po::variables_map &vm) override final
	Write from data to output file. More...

Protected Member Functions inherited from Nektar::FieldUtils::OutputVtkBase
virtual std::string	v_GetModuleName () override

virtual void	v_OutputFromPts (po::variables_map &vm) override
	Write from pts to output file. More...

virtual void	v_OutputFromExp (po::variables_map &vm) override
	Write from m_exp to output file. More...

virtual void	v_OutputFromData (po::variables_map &vm) override
	Write from data to output file. More...

virtual fs::path	v_GetPath (std::string &filename, po::variables_map &vm) override

virtual fs::path	v_GetFullOutName (std::string &filename, po::variables_map &vm) override

std::string	PrepareOutput (po::variables_map &vm)

Protected Member Functions inherited from Nektar::FieldUtils::OutputFileBase
virtual void	v_Process (po::variables_map &vm) override
	Write fld to output file. More...

virtual std::string	v_GetModuleName () override

virtual std::string	v_GetModuleDescription () override

virtual ModulePriority	v_GetModulePriority () override

virtual void	v_OutputFromPts (po::variables_map &vm)=0
	Write from pts to output file. More...

virtual void	v_OutputFromExp (po::variables_map &vm)=0
	Write from m_exp to output file. More...

virtual void	v_OutputFromData (po::variables_map &vm)=0
	Write from data to output file. More...

virtual fs::path	v_GetPath (std::string &filename, po::variables_map &vm)

fs::path	GetPath (std::string &filename, po::variables_map &vm)

virtual fs::path	v_GetFullOutName (std::string &filename, po::variables_map &vm)

fs::path	GetFullOutName (std::string &filename, po::variables_map &vm)

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 ()

Protected Attributes
vtkSmartPointer< vtkUnstructuredGrid >	m_vtkMesh
	Cache file for unstructured grid VTK mesh data. More...

int	m_numPlanes = 1
	Number of planes if homogeneous. More...

bool	m_extraPlane = false
	Flag if extra plane in case of fourier expansion in homogeneous dir. More...

bool	m_cachedMesh = false
	Flag if mesh has been cached. More...

Protected Attributes inherited from Nektar::FieldUtils::OutputFileBase
bool	m_requireEquiSpaced

Protected Attributes inherited from Nektar::FieldUtils::OutputModule
std::ofstream	m_fldFile
	Output stream. 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...

Private Member Functions
vtkSmartPointer< vtkUnstructuredGrid >	OutputFromExpHighOrder (po::variables_map &vm)
	Prepare high order Lagrange VTK output. More...

void	AddFieldDataToVTKHighOrder (po::variables_map &vm, std::string &filename, vtkSmartPointer< vtkUnstructuredGrid > &vtkMesh)
	Add field data to high order Lagrange VTK output. More...

vtkSmartPointer< vtkUnstructuredGrid >	OutputFromExpLowOrder ()
	Prepare low order VTK output. More...

void	AddFieldDataToVTKLowOrder (po::variables_map &vm, std::string &filename, vtkSmartPointer< vtkUnstructuredGrid > &vtkMesh)
	Add field data to low order VTK output. More...

void	OutputFromExpLowOrderMultiBlock (po::variables_map &vm, std::string &filename)
	Prepare low order multi-block VTK output & add field data. More...

void	WriteVTK (vtkDataObject *vtkMesh, std::string &filename, po::variables_map &vm)
	Write VTK file using. More...

void	WritePVtu (po::variables_map &vm)
	Write the parallel .pvtu file. More...

Additional Inherited Members
Public Attributes inherited from Nektar::FieldUtils::Module
FieldSharedPtr	m_f
	Field object. More...

Detailed Description

Converter from fld to vtk.

Definition at line 51 of file OutputVtk.h.

Constructor & Destructor Documentation

◆ OutputVtk()

Nektar::FieldUtils::OutputVtk::OutputVtk ( FieldSharedPtr f )

explicit

Definition at line 59 of file OutputVtk.cpp.

                                     : OutputVtkBase(std::move(f))
{
    m_requireEquiSpaced = true;
    m_config["legacy"] =
        ConfigOption(true, "0", "Output using legacy manual file writers");
    m_config["highorder"] = ConfigOption(
        true, "0", "Output using new high-order Lagrange elements");
    m_config["multiblock"] = ConfigOption(
        true, "0", "Output using multi-blocks to separate composites");
    m_config["uncompress"] = ConfigOption(true, "0", "Uncompress xml sections");
}

References Nektar::FieldUtils::Module::m_config, and Nektar::FieldUtils::OutputFileBase::m_requireEquiSpaced.

◆ ~OutputVtk()

Nektar::FieldUtils::OutputVtk::~OutputVtk ( )

finaldefault

Member Function Documentation

◆ AddFieldDataToVTKHighOrder()

void Nektar::FieldUtils::OutputVtk::AddFieldDataToVTKHighOrder	(	po::variables_map &	vm,
		std::string &	filename,
		vtkSmartPointer< vtkUnstructuredGrid > &	vtkMesh
	)

private

Add field data to high order Lagrange VTK output.

Definition at line 1431 of file OutputVtk.cpp.

{
    // Convert expansion to an equispaced points field
    if (m_cachedMesh)
    {
        auto equispaced =
            MemoryManager<ProcessEquiSpacedOutput>::AllocateSharedPtr(m_f);
        equispaced->Process(vm);
    }
    LibUtilities::PtsFieldSharedPtr fPts = m_f->m_fieldPts;
 
    int nPts = static_cast<int>(fPts->GetNpoints());
    int dim  = static_cast<int>(fPts->GetDim());
 
    Array<OneD, Array<OneD, NekDouble>> pts;
    fPts->GetPts(pts);
 
    // Insert field information
    for (int i = 0; i < fPts->GetNFields(); ++i)
    {
        vtkNew<vtkDoubleArray> fieldData;
        fieldData->SetArray(&pts[dim + i][0], nPts, 1);
        fieldData->SetName(&fPts->GetFieldName(i)[0]);
        vtkMesh->GetPointData()->AddArray(fieldData.GetPointer());
    }
 
    WriteVTK(vtkMesh, filename, vm);
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), m_cachedMesh, Nektar::FieldUtils::Module::m_f, and WriteVTK().

Referenced by v_OutputFromExp().

◆ AddFieldDataToVTKLowOrder()

void Nektar::FieldUtils::OutputVtk::AddFieldDataToVTKLowOrder	(	po::variables_map &	vm,
		std::string &	filename,
		vtkSmartPointer< vtkUnstructuredGrid > &	vtkMesh
	)

private

Add field data to low order VTK output.

Definition at line 979 of file OutputVtk.cpp.

{
    // Insert field information we iterate by element from first plane to last,
    // while the fld file iterates by plane from first element to last
    int nElmts         = static_cast<int>(m_f->m_exp[0]->GetNumElmts());
    int numPtsPerPlane = m_f->m_exp[0]->GetTotPoints() / m_numPlanes;
    for (int v = 0; v < m_f->m_variables.size(); ++v)
    {
        vtkNew<vtkDoubleArray> fieldData;
        for (int i = 0; i < nElmts; ++i)
        {
            int elmtOffset          = m_f->m_exp[v]->GetPhys_Offset(i);
            int nPtsPerElmtPerPlane = m_f->m_exp[v]->GetExp(i)->GetTotPoints();
 
            for (int j = 0; j < m_numPlanes; ++j)
            {
                int planeOffset = j * numPtsPerPlane;
                for (int k = 0; k < nPtsPerElmtPerPlane; ++k)
                {
                    fieldData->InsertNextValue(
                        m_f->m_exp[v]->GetPhys()[elmtOffset + planeOffset + k]);
                }
            }
 
            // if extra plane we copy the first plane values in to last plane
            if (m_extraPlane)
            {
                for (int k = 0; k < nPtsPerElmtPerPlane; ++k)
                {
                    fieldData->InsertNextValue(
                        m_f->m_exp[v]->GetPhys()[elmtOffset + k]);
                }
            }
        }
 
        fieldData->SetName(&m_f->m_variables[v][0]);
        vtkMesh->GetPointData()->AddArray(fieldData.GetPointer());
    }
 
    WriteVTK(vtkMesh, filename, vm);
}

References m_extraPlane, Nektar::FieldUtils::Module::m_f, m_numPlanes, and WriteVTK().

Referenced by v_OutputFromExp().

◆ create()

static std::shared_ptr< Module > Nektar::FieldUtils::OutputVtk::create ( const FieldSharedPtr & f )

inlinestatic

Creates an instance of this class.

Definition at line 55 of file OutputVtk.h.

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

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

◆ OutputFromExpHighOrder()

vtkSmartPointer< vtkUnstructuredGrid > Nektar::FieldUtils::OutputVtk::OutputFromExpHighOrder ( po::variables_map & vm )

private

Prepare high order Lagrange VTK output.

Definition at line 1302 of file OutputVtk.cpp.

{
    ASSERTL0(
        m_f->m_numHomogeneousDir == 0,
        "High order VTK is not implemented for homogeneous expansion types.")
 
    // Save shapetypes before convert to equispaced because that nukes the
    // explist
    std::vector<LibUtilities::ShapeType> sType;
    for (auto &i : *m_f->m_exp[0]->GetExp())
    {
        sType.emplace_back(i->GetGeom()->GetShapeType());
    }
 
    // Convert expansion to an equispaced points field
    auto equispaced =
        MemoryManager<ProcessEquiSpacedOutput>::AllocateSharedPtr(m_f);
    equispaced->Process(vm);
 
    // Insert points
    vtkNew<vtkPoints> vtkPoints;
    LibUtilities::PtsFieldSharedPtr fPts = m_f->m_fieldPts;
    vtkPoints->SetDataType(VTK_DOUBLE);
 
    Array<OneD, Array<OneD, NekDouble>> pts;
    fPts->GetPts(pts);
 
    int dim  = static_cast<int>(fPts->GetDim());
    int nPts = static_cast<int>(fPts->GetNpoints());
    switch (dim)
    {
        case 3:
            for (int i = 0; i < nPts; ++i)
            {
                vtkPoints->InsertNextPoint(pts[0][i], pts[1][i], pts[2][i]);
            }
            break;
        case 2:
            for (int i = 0; i < nPts; ++i)
            {
                vtkPoints->InsertNextPoint(pts[0][i], pts[1][i], 0.0);
            }
            break;
        case 1:
            for (int i = 0; i < nPts; ++i)
            {
                vtkPoints->InsertNextPoint(pts[0][i], 0.0, 0.0);
            }
            break;
        default:
            break;
    }
 
    // Mesh file to output
    vtkSmartPointer<vtkUnstructuredGrid> vtkMesh =
        vtkSmartPointer<vtkUnstructuredGrid>::New();
    vtkMesh->SetPoints(vtkPoints.GetPointer());
 
    // Cache ordering for shape type & npts so we aren't recreating mappings
    std::unordered_map<size_t, std::vector<long long>> mappingCache;
 
    // Get offset per element in a vector from ppe
    std::vector<int> ppe = m_f->m_fieldPts->GetPointsPerElement();
    Array<OneD, int> ppeOffset(ppe.size() + 1, 0.0);
    std::partial_sum(ppe.begin(), ppe.end(), &ppeOffset[1]);
 
    // Insert elements
    for (int i = 0; i < ppe.size(); ++i)
    {
        // Construct inverse of input reordering.
        // First try to find it in our mapping cache.
        auto oIt = mappingCache.find(key2(sType[i], ppe[i]));
        if (oIt == mappingCache.end())
        {
            std::vector<long long> p(ppe[i]);
            std::iota(p.begin(), p.end(), 0);
            switch (sType[i])
            {
                case LibUtilities::eQuadrilateral:
                    p = quadTensorNodeOrdering(p);
                    break;
                case LibUtilities::eTriangle:
                    p = triTensorNodeOrdering(p);
                    break;
                case LibUtilities::eTetrahedron:
                    p = tetTensorNodeOrdering(p);
                    break;
                case LibUtilities::eHexahedron:
                    p = hexTensorNodeOrdering(p);
                    break;
                case LibUtilities::ePrism:
                    p = prismTensorNodeOrdering(p);
                    break;
                default:
                    NEKERROR(ErrorUtil::efatal,
                             "High-order VTU output not set up for the " +
                                 static_cast<std::string>(
                                     LibUtilities::ShapeTypeMap[sType[i]]) +
                                 " shape type.");
                    break;
            }
 
            // Invert the ordering as this is for spectral -> recursive (VTU)
            std::vector<long long> inv(ppe[i]);
            for (int j = 0; j < ppe[i]; ++j)
            {
                inv[p[j]] = j;
            }
 
            oIt =
                mappingCache.insert(std::make_pair(key2(sType[i], ppe[i]), inv))
                    .first;
        }
 
        // Add offset to reordering
        std::vector<long long> p = oIt->second;
        for (long long &j : p)
        {
            j += ppeOffset[i];
        }
 
        vtkMesh->InsertNextCell(GetHighOrderVtkCellType(sType[i]), ppe[i],
                                &p[0]);
    }
 
    return vtkMesh;
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::ErrorUtil::efatal, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTetrahedron, Nektar::LibUtilities::eTriangle, Nektar::FieldUtils::Module::m_f, NEKERROR, CellMLToNektar.cellml_metadata::p, and Nektar::LibUtilities::ShapeTypeMap.

Referenced by v_OutputFromExp().

◆ OutputFromExpLowOrder()

vtkSmartPointer< vtkUnstructuredGrid > Nektar::FieldUtils::OutputVtk::OutputFromExpLowOrder ( )

private

Prepare low order VTK output.

Definition at line 856 of file OutputVtk.cpp.

{
    // Mesh file to output
    vtkSmartPointer<vtkUnstructuredGrid> vtkMesh =
        vtkSmartPointer<vtkUnstructuredGrid>::New();
 
    // Cache ordering so we aren't recreating mappings
    std::unordered_map<size_t, std::vector<long long>> mappingCache;
 
    // Choose cell types and num quad points based on dimension
    int nDim     = m_f->m_exp[0]->GetShapeDimension();
    int nHomoDir = m_f->m_numHomogeneousDir;
 
    auto type = (nDim + nHomoDir == 3)   ? VTK_HEXAHEDRON
                : (nDim + nHomoDir == 2) ? VTK_QUAD
                                         : VTK_LINE;
    int nQpts = (nDim + nHomoDir == 3) ? 8 : (nDim + nHomoDir == 2) ? 4 : 2;
 
    // Fill homogeneous info
    if (nHomoDir != 0)
    {
        ASSERTL0(nHomoDir == 1 &&
                     m_f->m_exp[0]->GetExpType() == MultiRegions::e3DH1D,
                 "Only regular expansions and the 3DH1D homogeneous expansion "
                 "are supported in the new VTK writer. Please use the 'legacy' "
                 "option for all other expansion types.")
 
        m_numPlanes = m_f->m_exp[0]->GetHomogeneousBasis()->GetNumModes();
        // Extra plane if fourier
        m_extraPlane = (m_f->m_exp[0]->GetHomogeneousBasis()->GetBasisType() ==
                            LibUtilities::eFourier &&
                        m_f->m_exp[0]->GetHomogeneousBasis()->GetPointsType() ==
                            LibUtilities::eFourierEvenlySpaced);
    }
 
    // Insert points
    vtkNew<vtkPoints> vtkPoints;
    vtkPoints->SetDataType(VTK_DOUBLE);
 
    int offset = 0;
    int nElmts = static_cast<int>(m_f->m_exp[0]->GetNumElmts());
    for (int i = 0; i < nElmts; ++i)
    {
        Array<OneD, int> nquad(3, 0);
        for (int j = 0; j < nDim; ++j)
        {
            nquad[j] = m_f->m_exp[0]->GetExp(i)->GetNumPoints(j);
        }
 
        // Add homogeneous direction num points
        if (nHomoDir != 0)
        {
            nquad[nDim] = m_numPlanes;
 
            // Add extra plane if fourier
            if (m_extraPlane)
            {
                nquad[nDim]++;
            }
        }
 
        // Get total points by multiplying all nquads together (accounts for
        // homo)
        int nPts = nquad[0];
        for (int j = 1; j < nDim + nHomoDir; ++j)
        {
            nPts *= nquad[j];
        }
 
        // Get all coordinates at quadrature points in this element
        Array<OneD, NekDouble> x(nPts, 0.0), y(nPts, 0.0), z(nPts, 0.0);
        m_f->m_exp[0]->GetCoords(i, x, y, z);
 
        // If add extra plane for fourier need to fill last plane coordinates
        if (m_extraPlane)
        {
            // Copy x & y to extra plane
            Array<OneD, NekDouble> tmp;
            Vmath::Vcopy(nquad[0] * nquad[1], x, 1,
                         tmp = x + (nquad[2] - 1) * nquad[0] * nquad[1], 1);
            Vmath::Vcopy(nquad[0] * nquad[1], y, 1,
                         tmp = y + (nquad[2] - 1) * nquad[0] * nquad[1], 1);
 
            // Fill z on extra plane
            NekDouble zDouble = z[nquad[0] * nquad[1] * (nquad[2] - 1) - 1] +
                                (z[nquad[0] * nquad[1]] - z[0]);
            Vmath::Fill(nquad[0] * nquad[1], zDouble,
                        tmp = z + (nquad[2] - 1) * nquad[0] * nquad[1], 1);
        }
 
        // Insert points in to vtk mesh
        for (int j = 0; j < nPts; ++j)
        {
            vtkPoints->InsertNextPoint(x[j], y[j], z[j]);
        }
 
        // Insert elements
        auto oIt = mappingCache.find(key3(nquad[0], nquad[1], nquad[2]));
        if (oIt == mappingCache.end())
        {
            auto p = lowOrderMapping(nDim + nHomoDir, nquad);
            oIt    = mappingCache
                      .insert(
                          std::make_pair(key3(nquad[0], nquad[1], nquad[2]), p))
                      .first;
        }
 
        auto p = oIt->second;
        std::for_each(p.begin(), p.end(),
                      [offset](long long &d) { d += offset; });
        for (int j = 0; j < p.size(); j += nQpts)
        {
            vtkMesh->InsertNextCell(type, nQpts, &p[j]);
        }
 
        offset += nPts;
    }
 
    vtkMesh->SetPoints(vtkPoints.GetPointer());
 
    return vtkMesh;
}

References ASSERTL0, Nektar::UnitTests::d(), Nektar::MultiRegions::e3DH1D, Nektar::LibUtilities::eFourier, Nektar::LibUtilities::eFourierEvenlySpaced, Vmath::Fill(), m_extraPlane, Nektar::FieldUtils::Module::m_f, m_numPlanes, CellMLToNektar.cellml_metadata::p, Vmath::Vcopy(), and Nektar::UnitTests::z().

Referenced by v_OutputFromExp().

◆ OutputFromExpLowOrderMultiBlock()

void Nektar::FieldUtils::OutputVtk::OutputFromExpLowOrderMultiBlock	(	po::variables_map &	vm,
		std::string &	filename
	)

private

Prepare low order multi-block VTK output & add field data.

Definition at line 1023 of file OutputVtk.cpp.

{
    ASSERTL0(
        m_f->m_numHomogeneousDir == 0,
        "Multi block VTK is not implemented for homogeneous expansion types.")
 
    ASSERTL0(m_f->m_comm->GetSpaceComm()->IsSerial(),
             "Multi block VTK is not implemented in parallel.")
 
    int dim = m_f->m_graph->GetMeshDimension();
 
    // Create mappings from geometry id to expansion ids
    std::array<std::map<int, std::pair<int, int>>, 4> geomIdToExpId;
    int nElmts = static_cast<int>(m_f->m_exp[0]->GetNumElmts());
    for (int i = 0; i < nElmts; ++i)
    {
        auto geom = m_f->m_exp[0]->GetExp(i)->GetGeom();
        geomIdToExpId[geom->GetShapeDim()][geom->GetGlobalID()] =
            std::make_pair(i, -1);
 
        for (int j = 0; j < geom->GetNumFaces(); ++j)
        {
            geomIdToExpId[2][geom->GetFid(j)] = std::make_pair(i, j);
        }
 
        for (int j = 0; j < geom->GetNumEdges(); ++j)
        {
            geomIdToExpId[1][geom->GetEid(j)] = std::make_pair(i, j);
        }
    }
 
    // Cache ordering so we aren't recreating mappings
    std::unordered_map<size_t, std::vector<long long>> mappingCache;
 
    std::map<int, vtkNew<vtkUnstructuredGrid>> vtkMesh;
    std::map<int, SpatialDomains::CompositeSharedPtr> composites =
        m_f->m_graph->GetComposites();
    std::map<int, std::string> compositeNames;
    for (auto &comp : composites)
    {
        // Vector of field data
        std::vector<vtkNew<vtkDoubleArray>> fieldData(m_f->m_variables.size());
 
        // Insert points
        vtkNew<vtkPoints> vtkPoints;
        vtkPoints->SetDataType(VTK_DOUBLE);
 
        int compId = comp.first;
        std::vector<std::shared_ptr<SpatialDomains::Geometry>> geomVec =
            comp.second->m_geomVec;
 
        unsigned int offset = 0;
        for (auto &geom : geomVec)
        {
            int geomId = geom->GetGlobalID();
            int sDim   = geom->GetShapeDim();
            auto type  = (sDim == 3)   ? VTK_HEXAHEDRON
                         : (sDim == 2) ? VTK_QUAD
                                       : VTK_LINE;
            int nQpts  = (sDim == 3) ? 8 : (sDim == 2) ? 4 : 2;
 
            LocalRegions::ExpansionSharedPtr exp =
                m_f->m_exp[0]->GetExp(geomIdToExpId[sDim][geomId].first);
 
            unsigned int nPts = exp->GetTotPoints();
            Array<OneD, NekDouble> x(nPts, 0.0), y(nPts, 0.0), z(nPts, 0.0);
            exp->GetCoords(x, y, z);
 
            int offsetPhys = m_f->m_exp[0]->GetPhys_Offset(
                geomIdToExpId[sDim][geomId].first);
            if (sDim == dim)
            {
                for (int j = 0; j < nPts; ++j)
                {
                    vtkPoints->InsertNextPoint(x[j], y[j], z[j]);
 
                    // Add field data
                    for (int k = 0; k < m_f->m_variables.size(); ++k)
                    {
                        fieldData[k]->InsertNextValue(
                            m_f->m_exp[k]->GetPhys()[j + offsetPhys]);
                    }
                }
            }
            else
            {
                Array<OneD, int> pointsMap;
                exp->GetTracePhysMap(geomIdToExpId[sDim][geomId].second,
                                     pointsMap);
                for (int j : pointsMap)
                {
                    vtkPoints->InsertNextPoint(x[j], y[j], z[j]);
 
                    // Add field data
                    for (int k = 0; k < m_f->m_variables.size(); ++k)
                    {
                        fieldData[k]->InsertNextValue(
                            m_f->m_exp[k]->GetPhys()[offsetPhys + j]);
                    }
                }
 
                exp  = exp->GetTraceExp(geomIdToExpId[sDim][geomId].second);
                nPts = pointsMap.size();
            }
 
            Array<OneD, int> nquad(3, 0);
            for (int j = 0; j < sDim; ++j)
            {
                nquad[j] = exp->GetNumPoints(j);
            }
 
            auto oIt = mappingCache.find(key3(nquad[0], nquad[1], nquad[2]));
            if (oIt == mappingCache.end())
            {
                auto p = lowOrderMapping(sDim, nquad);
                oIt    = mappingCache
                          .insert(std::make_pair(
                              key3(nquad[0], nquad[1], nquad[2]), p))
                          .first;
            }
 
            auto p = oIt->second;
            std::for_each(p.begin(), p.end(),
                          [offset](long long &d) { d += offset; });
 
            for (int j = 0; j < p.size(); j += nQpts)
            {
                vtkMesh[compId]->InsertNextCell(type, nQpts, &p[j]);
            }
 
            offset += nPts;
        }
 
        vtkMesh[compId]->SetPoints(vtkPoints.GetPointer());
 
        // Add all fields to vtkMesh
        for (int k = 0; k < m_f->m_variables.size(); ++k)
        {
            fieldData[k]->SetName(&m_f->m_variables[k][0]);
            vtkMesh[compId]->GetPointData()->AddArray(
                fieldData[k].GetPointer());
        }
 
        // Find composite names if exists and store in map
        // Set name as boundary label if it exists otherwise index ID
        compositeNames[compId] = "Composite ID " + std::to_string(compId);
        auto clabels           = m_f->m_graph->GetCompositesLabels();
        auto oIt               = clabels.find(compId);
        if (oIt != clabels.end())
        {
            compositeNames[compId] = oIt->second;
        }
    }
 
    // Main multi-block set
    vtkNew<vtkMultiBlockDataSet> mainBlock;
 
    std::set<int> compSet;
 
    // Fill domain multi blocks from composites
    vtkNew<vtkMultiBlockDataSet> mainDomainBlock;
    auto domains = m_f->m_graph->GetDomain();
    std::vector<vtkNew<vtkMultiBlockDataSet>> domainMultiBlocks(domains.size());
    for (int i = 0; i < domains.size(); ++i)
    {
        auto dom = domains[i];
 
        // Loop over composites and see if in domain
        for (auto &comp : composites)
        {
            int compId = comp.first;
            if (dom.find(compId) != dom.end())
            {
                unsigned int nBlock = domainMultiBlocks[i]->GetNumberOfBlocks();
                domainMultiBlocks[i]->SetBlock(nBlock,
                                               vtkMesh[compId].GetPointer());
                domainMultiBlocks[i]->GetMetaData(nBlock)->Set(
                    vtkCompositeDataSet::NAME(),
                    compositeNames[compId].c_str());
                compSet.insert(compId);
            }
        }
 
        unsigned int nBlock = mainDomainBlock->GetNumberOfBlocks();
        mainDomainBlock->SetBlock(nBlock, domainMultiBlocks[i].GetPointer());
        mainDomainBlock->GetMetaData(nBlock)->Set(
            vtkCompositeDataSet::NAME(),
            ("Domain ID " + std::to_string(i)).c_str());
    }
 
    if (mainDomainBlock->GetNumberOfBlocks() != 0)
    {
        auto nBlock = static_cast<unsigned int>(mainBlock->GetNumberOfBlocks());
        mainBlock->SetBlock(nBlock, mainDomainBlock.GetPointer());
        mainBlock->GetMetaData(nBlock)->Set(vtkCompositeDataSet::NAME(),
                                            "Domains");
    }
 
    // Fill boundary multi blocks from composites
    SpatialDomains::BoundaryConditions bcs(m_f->m_session,
                                           m_f->m_exp[0]->GetGraph());
    const SpatialDomains::BoundaryRegionCollection &bregions =
        bcs.GetBoundaryRegions();
 
    vtkNew<vtkMultiBlockDataSet> mainBoundaryBlock;
    std::vector<vtkNew<vtkMultiBlockDataSet>> boundaryMultiBlocks(
        bregions.size());
    int cnt = 0;
    for (auto &boundary : bregions)
    {
        // Loop over composites and see if in boundary
        for (auto &comp : composites)
        {
            int compId = comp.first;
            if (boundary.second->find(compId) != boundary.second->end())
            {
                unsigned int nBlock =
                    boundaryMultiBlocks[cnt]->GetNumberOfBlocks();
                boundaryMultiBlocks[cnt]->SetBlock(
                    nBlock, vtkMesh[compId].GetPointer());
                boundaryMultiBlocks[cnt]->GetMetaData(nBlock)->Set(
                    vtkCompositeDataSet::NAME(),
                    compositeNames[compId].c_str());
                compSet.insert(compId);
            }
        }
 
        unsigned int nBlock = mainBoundaryBlock->GetNumberOfBlocks();
        mainBoundaryBlock->SetBlock(nBlock,
                                    boundaryMultiBlocks[cnt++].GetPointer());
 
        // Set name as boundary label if it exists otherwise index ID
        std::string name = "Boundary ID " + std::to_string(boundary.first);
        auto blabels     = bcs.GetBoundaryLabels();
        auto oIt         = blabels.find(boundary.first);
        if (oIt != blabels.end())
        {
            name = oIt->second;
        }
 
        mainBoundaryBlock->GetMetaData(nBlock)->Set(vtkCompositeDataSet::NAME(),
                                                    name.c_str());
    }
 
    if (mainBoundaryBlock->GetNumberOfBlocks() != 0)
    {
        auto nBlock = static_cast<unsigned int>(mainBlock->GetNumberOfBlocks());
        mainBlock->SetBlock(nBlock, mainBoundaryBlock.GetPointer());
        mainBlock->GetMetaData(nBlock)->Set(vtkCompositeDataSet::NAME(),
                                            "Boundaries");
    }
 
    // Include all other composites (not domains or boundaries)
    vtkNew<vtkMultiBlockDataSet> compsMultiBlocks;
    cnt = 0;
    for (auto &comp : composites)
    {
        int compId = comp.first;
        if (compSet.find(compId) == compSet.end())
        {
            unsigned int nBlock = compsMultiBlocks->GetNumberOfBlocks();
            compsMultiBlocks->SetBlock(nBlock, vtkMesh[compId].GetPointer());
            compsMultiBlocks->GetMetaData(nBlock)->Set(
                vtkCompositeDataSet::NAME(), compositeNames[compId].c_str());
        }
    }
 
    if (compsMultiBlocks->GetNumberOfBlocks() != 0)
    {
        auto nBlock = static_cast<unsigned int>(mainBlock->GetNumberOfBlocks());
        mainBlock->SetBlock(nBlock, compsMultiBlocks.GetPointer());
        mainBlock->GetMetaData(nBlock)->Set(vtkCompositeDataSet::NAME(),
                                            "Other composites");
    }
 
    WriteVTK(mainBlock.GetPointer(), filename, vm);
}

References ASSERTL0, Nektar::UnitTests::d(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryLabels(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), Nektar::FieldUtils::Module::m_f, CellMLToNektar.pycml::name, CellMLToNektar.cellml_metadata::p, WriteVTK(), and Nektar::UnitTests::z().

Referenced by v_OutputFromExp().

◆ v_GetModuleName()

virtual std::string Nektar::FieldUtils::OutputVtk::v_GetModuleName ( )

inlinefinaloverrideprotectedvirtual

Reimplemented from Nektar::FieldUtils::OutputVtkBase.

Definition at line 67 of file OutputVtk.h.

    {
        return "OutputVtk";
    }

◆ v_OutputFromData()

void Nektar::FieldUtils::OutputVtk::v_OutputFromData ( po::variables_map & vm )

finaloverrideprotectedvirtual

Write from data to output file.

Reimplemented from Nektar::FieldUtils::OutputVtkBase.

Definition at line 1596 of file OutputVtk.cpp.

{
    boost::ignore_unused(vm);
    NEKERROR(ErrorUtil::efatal,
             "OutputVtk can't write using only FieldData. You may need "
             "to add a mesh XML file to your input files.");
}

References Nektar::ErrorUtil::efatal, and NEKERROR.

◆ v_OutputFromExp()

void Nektar::FieldUtils::OutputVtk::v_OutputFromExp ( po::variables_map & vm )

finaloverrideprotectedvirtual

Write from m_exp to output file.

Reimplemented from Nektar::FieldUtils::OutputVtkBase.

Definition at line 1609 of file OutputVtk.cpp.

{
    if (m_config["legacy"].m_beenSet)
    {
        ASSERTL0(!m_config["multiblock"].m_beenSet,
                 "Multi block VTK is not implemented for legacy output.")
 
        ASSERTL0(!m_config["highorder"].m_beenSet,
                 "High order VTK is not implemented for legacy output.")
 
        // No caching of mesh data in legacy output
        OutputVtkBase::v_OutputFromExp(vm);
        return;
    }
 
    // Extract the output filename and extension
    std::string filename = OutputVtkBase::PrepareOutput(vm);
 
    // Save mesh state (if using filter this allows us to only ProcessEquispaced
    // if needed)
    m_cachedMesh = m_vtkMesh;
 
    if (m_config["highorder"].m_beenSet)
    {
        ASSERTL0(!m_config["multiblock"].m_beenSet,
                 "Multi block VTK is not implemented for high-order output.")
 
        if (!m_cachedMesh)
        {
            m_vtkMesh = OutputFromExpHighOrder(vm);
        }
 
        AddFieldDataToVTKHighOrder(vm, filename, m_vtkMesh);
    }
    else if (m_config["multiblock"].m_beenSet)
    {
        // No mesh caching for multiblock due to complexity of field data
        OutputFromExpLowOrderMultiBlock(vm, filename);
    }
    else
    {
        if (!m_cachedMesh)
        {
            m_vtkMesh = OutputFromExpLowOrder();
        }
 
        AddFieldDataToVTKLowOrder(vm, filename, m_vtkMesh);
    }
}

References AddFieldDataToVTKHighOrder(), AddFieldDataToVTKLowOrder(), ASSERTL0, m_cachedMesh, Nektar::FieldUtils::Module::m_config, m_vtkMesh, OutputFromExpHighOrder(), OutputFromExpLowOrder(), OutputFromExpLowOrderMultiBlock(), Nektar::FieldUtils::OutputVtkBase::PrepareOutput(), and Nektar::FieldUtils::OutputVtkBase::v_OutputFromExp().

◆ v_OutputFromPts()

void Nektar::FieldUtils::OutputVtk::v_OutputFromPts ( po::variables_map & vm )

finaloverrideprotectedvirtual

Write from pts to output file.

Reimplemented from Nektar::FieldUtils::OutputVtkBase.

Definition at line 1604 of file OutputVtk.cpp.

{
    OutputVtkBase::v_OutputFromPts(vm);
}

References Nektar::FieldUtils::OutputVtkBase::v_OutputFromPts().

◆ WritePVtu()

void Nektar::FieldUtils::OutputVtk::WritePVtu ( po::variables_map & vm )

private

Write the parallel .pvtu file.

Definition at line 1524 of file OutputVtk.cpp.

{
    std::string filename = m_config["outfile"].as<std::string>();
    int dot              = static_cast<int>(filename.find_last_of('.'));
    std::string body     = filename.substr(0, dot);
    filename             = body + ".pvtu";
 
    std::ofstream outfile(filename.c_str());
 
    int nprocs = m_f->m_comm->GetSpaceComm()->GetSize();
    std::string path =
        LibUtilities::PortablePath(OutputVtkBase::v_GetPath(filename, vm));
 
    outfile << "<?xml version=\"1.0\"?>" << endl;
    outfile << "<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" "
            << "byte_order=\"LittleEndian\">" << endl;
    outfile << "  <PUnstructuredGrid GhostLevel=\"0\">" << endl;
 
    // Add time if exists
    if (!m_f->m_fieldMetaDataMap["Time"].empty())
    {
        outfile << "    <FieldData> " << endl;
        outfile << "      <DataArray type=\"Float64\" Name=\"TimeValue\" "
                   "NumberOfTuples=\"1\" format=\"ascii\">"
                << endl;
        outfile << m_f->m_fieldMetaDataMap["Time"] << std::endl;
        outfile << "      </DataArray>" << std::endl;
        outfile << "    </FieldData> " << endl;
    }
 
    // Add point coordinates
    outfile << "    <PPoints> " << endl;
    outfile << "      <PDataArray type=\"Float64\" NumberOfComponents=\"" << 3
            << "\"/> " << endl;
    outfile << "    </PPoints>" << endl;
 
    // Add cell information
    outfile << "    <PCells>" << endl;
    outfile << "      <PDataArray type=\"Int64\" Name=\"connectivity\" "
               "NumberOfComponents=\"1\"/>"
            << endl;
    outfile << "      <PDataArray type=\"Int64\" Name=\"offsets\"      "
               "NumberOfComponents=\"1\"/>"
            << endl;
    outfile << "      <PDataArray type=\"UInt8\" Name=\"types\"        "
               "NumberOfComponents=\"1\"/>"
            << endl;
    outfile << "    </PCells>" << endl;
 
    // Add fields (point data)
    outfile << "    <PPointData>" << endl;
    for (auto &var : m_f->m_variables)
    {
        outfile << "      <PDataArray type=\"Float64\" Name=\"" << var << "\"/>"
                << endl;
    }
    outfile << "    </PPointData>" << endl;
 
    // Add parallel files
    for (int i = 0; i < nprocs; ++i)
    {
        boost::format pad("P%1$07d.vtu");
        pad % i;
        outfile << "    <Piece Source=\"" << path << "/" << pad.str() << "\"/>"
                << endl;
    }
    outfile << "  </PUnstructuredGrid>" << endl;
    outfile << "</VTKFile>" << endl;
 
    cout << "Written file: " << filename << endl;
}

References CellMLToNektar.pycml::format, Nektar::FieldUtils::Module::m_config, Nektar::FieldUtils::Module::m_f, Nektar::LibUtilities::PortablePath(), and Nektar::FieldUtils::OutputVtkBase::v_GetPath().

Referenced by WriteVTK().

◆ WriteVTK()

void Nektar::FieldUtils::OutputVtk::WriteVTK	(	vtkDataObject *	vtkMesh,
		std::string &	filename,
		po::variables_map &	vm
	)

private

Write VTK file using.

Parameters

vtkMesh

Definition at line 1462 of file OutputVtk.cpp.

{
    // Initialise base writer class as default structured grid writer
    vtkSmartPointer<vtkXMLWriter> vtkMeshWriter =
        vtkNew<vtkXMLUnstructuredGridWriter>().GetPointer();
 
    // If multi block we need to switch to the MultiBlock writer
    if (m_config["multiblock"].m_beenSet)
    {
        vtkMeshWriter = vtkNew<vtkXMLMultiBlockDataWriter>().GetPointer();
    }
 
    // Choose the correct extension based on which mesh writer is used
    int dot          = static_cast<int>(filename.find_last_of('.'));
    std::string body = filename.substr(0, dot + 1);
    filename         = body + vtkMeshWriter->GetDefaultFileExtension();
 
    // Set time information if exists
    if (!m_f->m_fieldMetaDataMap["Time"].empty())
    {
        vtkMesh->GetInformation()->Set(
            vtkDataObject::DATA_TIME_STEP(),
            std::stod(m_f->m_fieldMetaDataMap["Time"]));
    }
    else if (!m_f->m_fieldMetaDataMap["FinalTime"].empty())
    {
        vtkMesh->GetInformation()->Set(
            vtkDataObject::DATA_TIME_STEP(),
            std::stod(m_f->m_fieldMetaDataMap["FinalTime"]));
    }
 
    // Write out the new mesh in XML format (don't support legacy
    // format here as we still have standard OutputVtk.cpp)
    vtkMeshWriter->SetFileName(filename.c_str());
 
#if VTK_MAJOR_VERSION <= 5
    vtkMeshWriter->SetInput(vtkMesh);
#else
    vtkMeshWriter->SetInputData(vtkMesh);
#endif
 
    if (m_config["uncompress"].m_beenSet)
    {
        vtkMeshWriter->SetDataModeToAscii();
    }
 
    vtkMeshWriter->Update();
 
    std::cout << "Written file: " << filename << std::endl;
 
    // Output parallel file information if needed - using a manual write.
    // We could use the VTK lib to do this, but that requires VTK with MPI
    // enabled & messing about with the parallel controller & changing
    // our file naming scheme as VTK forces _${proc-number} as a suffix...
    if (m_f->m_comm->GetSpaceComm()->TreatAsRankZero() &&
        !m_f->m_comm->GetSpaceComm()->IsSerial())
    {
        WritePVtu(vm);
    }
}