Extended ExportVTK to multigroup writing. Restructured Solveroutput for PN and...

Extended ExportVTK to multigroup writing. Restructured Solveroutput for PN and MN Solvers to Multigroup with cleaner code. Removed dFlux as Screenoutput. Adapted other solvers to maintain compatibility with the new ExportVTK format. _solverOutput is now legacy and will be removed in future commits!


Former-commit-id: 7b3475dd

code/include/solvers/mnsolver.h

@@ -24,8 +24,8 @@ class MNSolver : public Solver
     void Save( int currEnergy ) const override; /*! @brief Save Output solution at given energy (pseudo time) to VTK file */
 
   private:
-    unsigned _nTotalEntries;          /*! @brief: Total number of equations in the system */
-    unsigned short _nMaxMomentsOrder; /*! @brief: Max Order of Moments */
+    unsigned _nTotalEntries;    /*! @brief: Total number of equations in the system */
+    unsigned short _LMaxDegree; /*! @brief: Max Order of Moments */
 
     // Solver specific physics
     VectorVector _sigmaA; /*!  @brief: Absorption coefficient for all energies*/
@@ -50,7 +50,14 @@ class MNSolver : public Solver
     OptimizerBase* _optimizer; /*! @brief: Class to solve minimal entropy problem */
 
     // Output related members
-    std::vector<std::vector<double>> _outputFields; /*! @brief: Protoype output for multiple output fields. Will replace _solverOutput */
+    std::vector<std::vector<std::vector<double>>> _outputFields; /*! @brief: Solver Output: dimensions (GroupID,FieldID,CellID). !Protoype output for
+                                                                    multiple output fields. Will replace _solverOutput */
+    std::vector<std::vector<std::string>> _outputFieldNames;     /*! @brief: Names of the outputFields: dimensions (GroupID,FieldID) */
+
+    // ---- Member functions ---
+
+    /*! @brief Initializes the output groups and fields of this solver and names the fields */
+    void PrepareOutputFields();
 
     /*! @brief Function that writes NN Training Data in a .csv file */
     void WriteNNTrainingData( unsigned idx_pseudoTime );
@@ -60,7 +67,6 @@ class MNSolver : public Solver
     */
     double WriteOutputFields();
 
-    // Member functions
     /*! @brief : computes the global index of the moment corresponding to basis function (l,k)
      *  @param : degree l, it must hold: 0 <= l <=_nq
      *  @param : order k, it must hold: -l <=k <= l

code/include/solvers/pnsolver.h

@@ -42,8 +42,22 @@ class PNSolver : public Solver
     Vector _scatterMatDiag; /*! @brief: diagonal of the scattering matrix (its a diagonal matrix by construction) */
 
     // Output related members
-    std::vector<std::vector<double>> _outputFields; /*! @brief: Protoype output for multiple output fields. Will replace _solverOutput */
+    std::vector<std::vector<std::vector<double>>> _outputFields; /*! @brief: Solver Output: dimensions (GroupID,FieldID,CellID). !Protoype output for
+                                                                    multiple output fields. Will replace _solverOutput */
+    std::vector<std::vector<std::string>> _outputFieldNames;     /*! @brief: Names of the outputFields: dimensions (GroupID,FieldID) */
 
+    // ---- Member functions ----
+
+    // IO
+    /*! @brief Initializes the output groups and fields of this solver and names the fields */
+    void PrepareOutputFields();
+
+    /*! @brief Function that prepares VTK export and csv export of the current solver iteration
+        @returns: Mass of current iteration
+    */
+    double WriteOutputFields();
+
+    // Solver
     /*! @brief: parameter functions for setting up system matrix
      *  @param: degree l, it must hold: 0 <= l <=_nq
      *  @param : order k, it must hold: -l <=k <= l

code/include/solvers/solverbase.h

@@ -50,7 +50,7 @@ class Solver
 
     // Solution related members
     VectorVector _sol;                 /*! @brief solution of the PDE, e.g. angular flux or moments */
-    std::vector<double> _solverOutput; /*! @brief PROTOTYPE: Outputfield for solver ==> Will be replaced by _outputFields in the near future */
+    std::vector<double> _solverOutput; /*! @brief LEGACY: Outputfield for solver ==> Will be replaced by _outputFields in the near future */
 
     // we will have to add a further dimension for quadPoints and weights once we start with multilevel SN
 

code/src/common/io.cpp

@@ -37,8 +37,8 @@ using vtkCellDataToPointDataSP    = vtkSmartPointer<vtkCellDataToPointData>;
 // using vtkPointDataToCellDataSP    = vtkSmartPointer<vtkPointDataToCellData>;
 
 void ExportVTK( const std::string fileName,
-                const std::vector<std::vector<std::vector<double>>>& results,
-                const std::vector<std::string> fieldNames,
+                const std::vector<std::vector<std::vector<double>>>& outputFields,
+                const std::vector<std::vector<std::string>>& outputFieldNames,
                 const Mesh* mesh ) {
     int rank;
     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
@@ -96,39 +96,18 @@ void ExportVTK( const std::string fileName,
             grid->SetCells( VTK_QUAD, cellArray );
         }
 
-        for( unsigned idx_group = 0; idx_group < results.size(); idx_group++ ) {
-            // why is there a group dimension for results, but not for names?
-            // Right now, this is ok, since we only have one group... TODO!
+        // Write the output
+        for( unsigned idx_group = 0; idx_group < outputFields.size(); idx_group++ ) {
 
-            for( unsigned idx_field = 0; idx_field < results[idx_group].size(); idx_field++ ) {    // Loop over all output fields
+            for( unsigned idx_field = 0; idx_field < outputFields[idx_group].size(); idx_field++ ) {    // Loop over all output fields
 
                 auto cellData = vtkDoubleArraySP::New();
-                cellData->SetName( fieldNames[idx_field].c_str() );
+                cellData->SetName( outputFieldNames[idx_group][idx_field].c_str() );
 
                 for( unsigned idx_cell = 0; idx_cell < numCells; idx_cell++ ) {
-                    cellData->InsertNextValue( results[idx_group][idx_field][idx_cell] );
+                    cellData->InsertNextValue( outputFields[idx_group][idx_field][idx_cell] );
                 }
 
-                /*
-                switch( results[idx_group].size() ) {
-                    case 1:
-                        for( unsigned j = 0; j < numCells; j++ ) {
-                            cellData->InsertNextValue( results[idx_group][0][j] );
-                        }
-                        break;
-                    default:
-                        // for( unsigned l = 0; l < results[i].size(); l++ ) {    // Loop over all output fields
-                        //    for( unsigned j = 0; j < numCells; j++ ) {
-                        //        cellData->InsertNextValue( results[i][l][j] );
-                        //    }
-                        // }
-
-                        auto log = spdlog::get( "event" );
-                        ErrorMessages::Error( "Please implement output for results of size " + std::to_string( results[idx_group].size() ) + "!",
-                                              CURRENT_FUNCTION );
-                        break;
-                }*/
-
                 grid->GetCellData()->AddArray( cellData );
             }
         }

code/src/main.cpp

@@ -41,7 +41,6 @@ int main( int argc, char** argv ) {
     solver->Solve();
     solver->Save();
 
-    //  if( Py_IsInitialized() ) Py_Finalize();
     MPI_Finalize();
     return EXIT_SUCCESS;
 }

code/src/solvers/csdsnsolver.cpp

@@ -147,14 +147,16 @@ void CSDSNSolver::Solve() {
 
 void CSDSNSolver::Save() const {
     std::vector<std::string> fieldNames{ "dose" };
+    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };
     std::vector<std::vector<double>> scalarField( 1, _dose );
     std::vector<std::vector<std::vector<double>>> results{ scalarField };
-    ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
+    ExportVTK( _settings->GetOutputFile(), results, fieldNamesWrapper, _mesh );
 }
 
 void CSDSNSolver::Save( int currEnergy ) const {
-    std::vector<std::string> fieldNames{ "flux" };
+    std::vector<std::string> fieldNames{ "dose" };
+    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };
     std::vector<std::vector<double>> scalarField( 1, _solverOutput );
     std::vector<std::vector<std::vector<double>>> results{ scalarField };
-    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNames, _mesh );
+    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNamesWrapper, _mesh );
 }

code/src/solvers/mnsolver.cpp

@@ -6,6 +6,7 @@
 #include "optimizers/optimizerbase.h"
 #include "quadratures/quadraturebase.h"
 #include "solvers/sphericalharmonics.h"
+#include "toolboxes/errormessages.h"
 #include "toolboxes/textprocessingtoolbox.h"
 
 // externals
@@ -18,8 +19,8 @@
 MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
 
     // Is this good (fast) code using a constructor list?
-    _nMaxMomentsOrder = settings->GetMaxMomentDegree();
-    _nTotalEntries    = GlobalIndex( _nMaxMomentsOrder, int( _nMaxMomentsOrder ) ) + 1;
+    _LMaxDegree    = settings->GetMaxMomentDegree();
+    _nTotalEntries = GlobalIndex( _LMaxDegree, int( _LMaxDegree ) ) + 1;
 
     // build quadrature object and store quadrature points and weights
     _quadPoints       = _quadrature->GetPoints();
@@ -52,13 +53,13 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
     _alpha = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
 
     // Initialize and Pre-Compute Moments at quadrature points
-    _basis = new SphericalHarmonics( _nMaxMomentsOrder );
+    _basis = new SphericalHarmonics( _LMaxDegree );
 
     _moments = VectorVector( _nq, Vector( _nTotalEntries, 0.0 ) );
     ComputeMoments();
 
     // Solver output
-    _outputFields = std::vector( _nTotalEntries, std::vector( _nCells, 0.0 ) );
+    PrepareOutputFields();
 }
 
 MNSolver::~MNSolver() {
@@ -188,30 +189,79 @@ void MNSolver::Solve() {
     }
 }
 
+void MNSolver::PrepareOutputFields() {
+    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
+
+    _outputFieldNames.resize( nGroups );
+    _outputFields.resize( nGroups );
+
+    // Prepare all OutputGroups ==> Specified in option VOLUME_OUTPUT
+    for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
+        // Prepare all Output Fields per group
+
+        // Different procedure, depending on the Group...
+        switch( _settings->GetVolumeOutput()[idx_group] ) {
+            case MINIMAL:
+                // Currently only one entry ==> rad flux
+                _outputFields[idx_group].resize( 1 );
+                _outputFieldNames[idx_group].resize( 1 );
+
+                _outputFields[0][0].resize( _nCells );
+                _outputFieldNames[0][0] = "radiation flux density";
+                break;
+
+            case MOMENTS:
+                // As many entries as there are moments in the system
+                _outputFields[idx_group].resize( _nTotalEntries );
+                _outputFieldNames[idx_group].resize( _nTotalEntries );
+
+                for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
+                    for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
+                        _outputFields[idx_group][GlobalIndex( idx_l, idx_k )].resize( _nCells );
+
+                        _outputFieldNames[idx_group][GlobalIndex( idx_l, idx_k )] =
+                            std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
+                    }
+                }
+                break;
+            default: ErrorMessages::Error( "Volume Output Group not defined for MN Solver!", CURRENT_FUNCTION ); break;
+        }
+    }
+}
+
 double MNSolver::WriteOutputFields() {
-    double mass = 0.0;
+    double mass      = 0.0;
+    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
 
-    for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
-        for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
-            mass += _sol[idx_cell][0] * _areas[idx_cell];
-            _outputFields[idx_sys][idx_cell] = _sol[idx_cell][idx_sys];
-        }
+    // Compute total "mass" of the system ==> to check conservation properties
+    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
+        mass += _sol[idx_cell][0] * _areas[idx_cell];    // Should probably go to postprocessing
     }
 
+    for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
+        switch( _settings->GetVolumeOutput()[idx_group] ) {
+            case MINIMAL:
+                for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
+                    _outputFields[idx_group][0][idx_cell] = _sol[idx_cell][0];
+                }
+                break;
+            case MOMENTS:
+                for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
+                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
+                        _outputFields[idx_group][idx_sys][idx_cell] = _sol[idx_cell][idx_sys];
+                    }
+                }
+                break;
+            default: ErrorMessages::Error( "Volume Output Group not defined for MN Solver!", CURRENT_FUNCTION ); break;
+        }
+    }
     return mass;
 }
 
-void MNSolver::Save() const { Save( _nEnergies - 1 ); }
+void MNSolver::Save() const { ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); }
 
 void MNSolver::Save( int currEnergy ) const {
-    std::vector<std::string> fieldNames( _nTotalEntries, "flux" );
-    for( int idx_l = 0; idx_l <= (int)_nMaxMomentsOrder; idx_l++ ) {
-        for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
-            fieldNames[GlobalIndex( idx_l, idx_k )] = std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
-        }
-    }
-    std::vector<std::vector<std::vector<double>>> results{ _outputFields };
-    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNames, _mesh );
+    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), _outputFields, _outputFieldNames, _mesh );
 }
 
 void MNSolver::WriteNNTrainingData( unsigned idx_pseudoTime ) {

code/src/solvers/pnsolver.cpp

@@ -73,7 +73,7 @@ PNSolver::PNSolver( Config* settings ) : Solver( settings ) {
     // std::cout << "_nCells: " << _nCells << "\n";
 
     // Solver output
-    _outputFields = std::vector( _nTotalEntries, std::vector( _nCells, 0.0 ) );    // WILL REPLACE _solverOutput
+    PrepareOutputFields();
 }
 
 void PNSolver::Solve() {
@@ -83,30 +83,19 @@ void PNSolver::Solve() {
 
     auto log = spdlog::get( "event" );
 
-    // angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
     VectorVector psiNew = _sol;
-    double dFlux        = 1e10;
-    Vector fluxNew( _nCells, 0.0 );
-    Vector fluxOld( _nCells, 0.0 );
-
-    double mass1 = 0;
-    for( unsigned i = 0; i < _nCells; ++i ) {
-        _solverOutput[i] = _sol[i][0];
-        mass1 += _sol[i][0];
-    }
 
-    dFlux   = blaze::l2Norm( fluxNew - fluxOld );
-    fluxOld = fluxNew;
+    double mass = 0;
 
     Save( -1 );    // Save initial condition
 
     unsigned idx_system = 0;
 
-    if( rank == 0 ) log->info( "{:10}   {:10}", "t", "dFlux" );
-    if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", -1.0, dFlux, mass1 );
+    if( rank == 0 ) log->info( "{:10}   {:10}", "t", "mass" );
 
-    // Loop over energies (pseudo-time of continuous slowing down approach)
+    // if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", -1.0, dFlux, mass1 );
 
+    // Loop over energies (pseudo-time of continuous slowing down approach)
     for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
         // Loop over all spatial cells
         for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
@@ -152,27 +141,16 @@ void PNSolver::Solve() {
                 }
             }
         }
+
         // Update Solution
         _sol = psiNew;
 
-        // pseudo time iteration output
-        double mass = 0.0;
-        for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
-            for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
-
-                fluxNew[idx_cell] = _sol[idx_cell][0];    // zeroth moment is raditation densitiy we are interested in
-
-                _solverOutput[idx_cell] = _sol[idx_cell][0];
-                mass += _sol[idx_cell][0] * _areas[idx_cell];
-                _outputFields[idx_sys][idx_cell] = _sol[idx_cell][idx_sys];
-            }
-        }
-
-        // Screen Output
-        dFlux   = blaze::l2Norm( fluxNew - fluxOld );
-        fluxOld = fluxNew;
-        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", _energies[idx_energy], dFlux, mass );
+        // --- VTK and CSV Output ---
+        mass = WriteOutputFields();
         Save( idx_energy );
+
+        // --- Screen Output ---
+        if( rank == 0 ) log->info( "{:03.8f}  {:01.5e}", _energies[idx_energy], mass );
     }
 }
 
@@ -376,17 +354,79 @@ double PNSolver::LegendrePoly( double x, int l ) {    // Legacy. TO BE DELETED
     }
 }
 
-void PNSolver::Save() const { Save( _nEnergies - 1 ); }
+void PNSolver::PrepareOutputFields() {
+    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
 
-void PNSolver::Save( int currEnergy ) const {
-    std::vector<std::string> fieldNames( _nTotalEntries, "flux" );
-    for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
-        for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
-            fieldNames[GlobalIndex( idx_l, idx_k )] = std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
+    _outputFieldNames.resize( nGroups );
+    _outputFields.resize( nGroups );
+
+    // Prepare all OutputGroups ==> Specified in option VOLUME_OUTPUT
+    for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
+        // Prepare all Output Fields per group
+
+        // Different procedure, depending on the Group...
+        switch( _settings->GetVolumeOutput()[idx_group] ) {
+            case MINIMAL:
+                // Currently only one entry ==> rad flux
+                _outputFields[idx_group].resize( 1 );
+                _outputFieldNames[idx_group].resize( 1 );
+
+                _outputFields[0][0].resize( _nCells );
+                _outputFieldNames[0][0] = "radiation flux density";
+                break;
+
+            case MOMENTS:
+                // As many entries as there are moments in the system
+                _outputFields[idx_group].resize( _nTotalEntries );
+                _outputFieldNames[idx_group].resize( _nTotalEntries );
+
+                for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
+                    for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
+                        _outputFields[idx_group][GlobalIndex( idx_l, idx_k )].resize( _nCells );
+
+                        _outputFieldNames[idx_group][GlobalIndex( idx_l, idx_k )] =
+                            std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
+                    }
+                }
+                break;
+            default: ErrorMessages::Error( "Volume Output Group not defined for PN Solver!", CURRENT_FUNCTION ); break;
+        }
+    }
+}
+
+double PNSolver::WriteOutputFields() {
+    double mass      = 0.0;
+    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
+
+    // Compute total "mass" of the system ==> to check conservation properties
+    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
+        mass += _sol[idx_cell][0] * _areas[idx_cell];    // Should probably go to postprocessing
+    }
+
+    for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
+        switch( _settings->GetVolumeOutput()[idx_group] ) {
+            case MINIMAL:
+                for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
+                    _outputFields[idx_group][0][idx_cell] = _sol[idx_cell][0];
+                }
+                break;
+            case MOMENTS:
+                for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
+                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
+                        _outputFields[idx_group][idx_sys][idx_cell] = _sol[idx_cell][idx_sys];
+                    }
+                }
+                break;
+            default: ErrorMessages::Error( "Volume Output Group not defined for PN Solver!", CURRENT_FUNCTION ); break;
         }
     }
-    std::vector<std::vector<std::vector<double>>> results{ _outputFields };
-    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNames, _mesh );
+    return mass;
+}
+
+void PNSolver::Save() const { ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); }
+
+void PNSolver::Save( int currEnergy ) const {
+    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), _outputFields, _outputFieldNames, _mesh );
 }
 
 void PNSolver::CleanFluxMatrices() {

code/src/solvers/snsolver.cpp

@@ -153,18 +153,22 @@ void SNSolver::Solve() {
 
 void SNSolver::Save() const {
     std::vector<std::string> fieldNames{ "flux" };
+    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };
+
     std::vector<double> flux( _nCells, 0.0 );
     for( unsigned i = 0; i < _nCells; ++i ) {
         flux[i] = dot( _sol[i], _weights );
     }
     std::vector<std::vector<double>> scalarField( 1, flux );
     std::vector<std::vector<std::vector<double>>> results{ scalarField };
-    ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
+    ExportVTK( _settings->GetOutputFile(), results, fieldNamesWrapper, _mesh );
 }
 
 void SNSolver::Save( int currEnergy ) const {
     std::vector<std::string> fieldNames{ "flux" };
+    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };
+
     std::vector<std::vector<double>> scalarField( 1, _solverOutput );
     std::vector<std::vector<std::vector<double>>> results{ scalarField };
-    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNames, _mesh );
+    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNamesWrapper, _mesh );
 }

code/src/solvers/snsolverMPI.cpp

@@ -69,13 +69,15 @@ void SNSolverMPI::Solve() {
 
 void SNSolverMPI::Save() const {
     std::vector<std::string> fieldNames{ "flux" };
+    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };
+
     std::vector<double> flux( _nCells, 0.0 );
     for( unsigned i = 0; i < _nCells; ++i ) {
         flux[i] = dot( _sol[i], _weights );
     }
     std::vector<std::vector<double>> scalarField( 1, flux );
     std::vector<std::vector<std::vector<double>>> results{ scalarField };
-    ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
+    ExportVTK( _settings->GetOutputFile(), results, fieldNamesWrapper, _mesh );
     auto log = spdlog::get( "event" );
     log->info( "Result successfully exported to '{0}'!", _settings->GetOutputFile() );
 }

code/src/solvers/solverbase.cpp

@@ -80,6 +80,8 @@ Solver* Solver::Create( Config* settings ) {
 
 void Solver::Save() const {
     std::vector<std::string> fieldNames{ "flux" };
+    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };
+
     std::vector<double> flux;
     flux.resize( _nCells );
 
@@ -88,5 +90,5 @@ void Solver::Save() const {
     }
     std::vector<std::vector<double>> scalarField( 1, flux );
     std::vector<std::vector<std::vector<double>>> results{ scalarField };
-    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( _nEnergies ), results, fieldNames, _mesh );
+    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( _nEnergies ), results, fieldNamesWrapper, _mesh );
 }

code/tests/test_cases.cpp

@@ -78,6 +78,7 @@ TEST_CASE( "linesource_PN", "[validation_tests]" ) {
     auto reference = readVTKFile( "../tests/input/linesource_PN_reference.vtk" );
 
     double eps = 1e-3;
+
     REQUIRE( test.size() == reference.size() );
     for( unsigned i = 0; i < test.size(); ++i ) {
         REQUIRE( std::fabs( test[i] - reference[i] ) < eps );
@@ -98,6 +99,7 @@ TEST_CASE( "linesource_MN", "[validation_tests]" ) {
         auto reference = readVTKFile( "../tests/input/linesource_MN_Quad_reference.vtk" );
 
         double eps = 1e-3;
+
         REQUIRE( test.size() == reference.size() );
         for( unsigned i = 0; i < test.size(); ++i ) {
             REQUIRE( std::fabs( test[i] - reference[i] ) < eps );

code/tests/test_image_conversion.cpp

@@ -79,11 +79,13 @@ TEST_CASE( "convert image data to grayscale matrix", "[image I/O]" ) {
         }
 
         std::vector<std::string> fieldNames{ "CT Data" };
+        std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };
+
         std::vector<std::vector<double>> scalarField( 1, result );
         std::vector<std::vector<std::vector<double>>> results{ scalarField };
         std::string outputFile = config->GetOutputFile();
         if( !TextProcessingToolbox::StringEndsWith( outputFile, ".vtk" ) ) outputFile.append( ".vtk" );
-        ExportVTK( outputFile, results, fieldNames, mesh );
+        ExportVTK( outputFile, results, fieldNamesWrapper, mesh );
 
         REQUIRE( std::filesystem::exists( outputFile ) );