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

virtual ModulePriority	v_GetModulePriority () override

fs::path	GetPath (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 ()

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::MultiRegions::e3DH1D, Nektar::LibUtilities::eFourier, Nektar::LibUtilities::eFourierEvenlySpaced, Vmath::Fill(), m_extraPlane, Nektar::FieldUtils::Module::m_f, m_numPlanes, CellMLToNektar.cellml_metadata::p, and Vmath::Vcopy().

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::SpatialDomains::BoundaryConditions::GetBoundaryLabels(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), Nektar::FieldUtils::Module::m_f, CellMLToNektar.pycml::name, CellMLToNektar.cellml_metadata::p, and WriteVTK().

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