diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 3225154fd181..ab6c1a2f91f6 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -3642,6 +3642,19 @@ class CConfig {
     return (Kind_Solver == HEAT_EQUATION) || (Kind_Solver == DISC_ADJ_HEAT);
   }
 
+  /*!
+   * \brief Return true if a high order FEM solver is in use.
+   */
+  bool GetFEMSolver(void) const {
+    switch (Kind_Solver) {
+      case FEM_EULER: case FEM_NAVIER_STOKES: case FEM_RANS: case FEM_LES:
+      case DISC_ADJ_FEM_EULER: case DISC_ADJ_FEM_NS: case DISC_ADJ_FEM_RANS:
+        return true;
+      default:
+        return false;
+    }
+  }
+
   /*!
    * \brief Kind of Multizone Solver.
    * \return Governing equation that we are solving.
diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp
index f71afcaa67ab..ed149c7fda36 100644
--- a/Common/src/geometry/CPhysicalGeometry.cpp
+++ b/Common/src/geometry/CPhysicalGeometry.cpp
@@ -166,13 +166,7 @@ CPhysicalGeometry::CPhysicalGeometry(CConfig *config, unsigned short val_iZone,
 
   /*--- Determine whether or not a FEM discretization is used ---*/
 
-  const bool fem_solver = ((config->GetKind_Solver() == FEM_EULER)          ||
-                           (config->GetKind_Solver() == FEM_NAVIER_STOKES)  ||
-                           (config->GetKind_Solver() == FEM_RANS)           ||
-                           (config->GetKind_Solver() == FEM_LES)            ||
-                           (config->GetKind_Solver() == DISC_ADJ_FEM_EULER) ||
-                           (config->GetKind_Solver() == DISC_ADJ_FEM_NS)    ||
-                           (config->GetKind_Solver() == DISC_ADJ_FEM_RANS));
+  const bool fem_solver = config->GetFEMSolver();
 
   /*--- Initialize counters for local/global points & elements ---*/
 
@@ -11756,4 +11750,4 @@ void CPhysicalGeometry::SetWallDistance(const CConfig *config, CADTElemClass *Wa
         node[iPoint]->SetWall_Distance(dist);
     }
   }
-}
\ No newline at end of file
+}
diff --git a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp
index 2923ca4d9eb5..18d7da37d96d 100644
--- a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp
+++ b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp
@@ -115,6 +115,9 @@ CDiscAdjSinglezoneDriver::CDiscAdjSinglezoneDriver(char* confFile,
 
 CDiscAdjSinglezoneDriver::~CDiscAdjSinglezoneDriver(void) {
 
+  delete direct_iteration;
+  delete direct_output;
+
 }
 
 void CDiscAdjSinglezoneDriver::Preprocess(unsigned long TimeIter) {
diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp
index 02bf4d3f15de..82f39eac864a 100644
--- a/SU2_CFD/src/drivers/CDriver.cpp
+++ b/SU2_CFD/src/drivers/CDriver.cpp
@@ -649,16 +649,10 @@ void CDriver::Input_Preprocessing(CConfig **&config, CConfig *&driver_config) {
     }
   }
 
-  /*--- Determine whether or not the FEM solver is used, which decides the
- type of geometry classes that are instantiated. Only adapted for single-zone problems ---*/
-
-  fem_solver = ((config_container[ZONE_0]->GetKind_Solver() == FEM_EULER)          ||
-                (config_container[ZONE_0]->GetKind_Solver() == FEM_NAVIER_STOKES)  ||
-                (config_container[ZONE_0]->GetKind_Solver() == FEM_RANS)           ||
-                (config_container[ZONE_0]->GetKind_Solver() == FEM_LES)            ||
-                (config_container[ZONE_0]->GetKind_Solver() == DISC_ADJ_FEM_EULER) ||
-                (config_container[ZONE_0]->GetKind_Solver() == DISC_ADJ_FEM_NS)    ||
-                (config_container[ZONE_0]->GetKind_Solver() == DISC_ADJ_FEM_RANS));
+  /*--- Determine whether or not the FEM solver is used, which decides the type of
+   *    geometry classes that are instantiated. Only adapted for single-zone problems ---*/
+
+  fem_solver = config_container[ZONE_0]->GetFEMSolver();
 
   fsi = config_container[ZONE_0]->GetFSI_Simulation();
 }
diff --git a/SU2_CFD/src/output/output_structure_legacy.cpp b/SU2_CFD/src/output/output_structure_legacy.cpp
index c3427fb4e23c..e4f8fd4c745a 100644
--- a/SU2_CFD/src/output/output_structure_legacy.cpp
+++ b/SU2_CFD/src/output/output_structure_legacy.cpp
@@ -4162,10 +4162,7 @@ void COutputLegacy::SetRestart(CConfig *config, CGeometry *geometry, CSolver **s
 
   /*--- Determine whether or not the FEM solver is used, which decides the
    type of geometry classes that are instantiated. ---*/
-  bool fem_solver = ((config->GetKind_Solver() == FEM_EULER)         ||
-                (config->GetKind_Solver() == FEM_NAVIER_STOKES) ||
-                (config->GetKind_Solver() == FEM_RANS)          ||
-                (config->GetKind_Solver() == FEM_LES));
+  const bool fem_solver = config->GetFEMSolver();
 
   unsigned long nPointTotal = 0;
   if ( fem_solver ) {
@@ -11880,11 +11877,7 @@ void COutputLegacy::SetResult_Files_Parallel(CSolver *****solver_container,
   for (iZone = 0; iZone < val_nZone; iZone++) {
 
     /*--- Bool to distinguish between the FVM and FEM solvers. ---*/
-    unsigned short KindSolver = config[iZone]->GetKind_Solver();
-    bool fem_solver = ((KindSolver == FEM_EULER) ||
-                       (KindSolver == FEM_NAVIER_STOKES) ||
-                       (KindSolver == FEM_RANS) ||
-                       (KindSolver == FEM_LES));
+    const bool fem_solver = config[iZone]->GetFEMSolver();
 
     /*--- Get the file output format ---*/
 
@@ -20393,11 +20386,7 @@ void COutputLegacy::PrepareOffsets(CConfig *config, CGeometry *geometry) {
 
   /*--- Bool to distinguish between the FVM and FEM solvers. ---*/
 
-  unsigned short KindSolver = config->GetKind_Solver();
-  bool fem_solver = ((KindSolver == FEM_EULER) ||
-                     (KindSolver == FEM_NAVIER_STOKES) ||
-                     (KindSolver == FEM_RANS) ||
-                     (KindSolver == FEM_LES));
+  const bool fem_solver = config->GetFEMSolver();
 
   /*--- Reset point sorting counters ---*/
 
diff --git a/SU2_DOT/src/SU2_DOT.cpp b/SU2_DOT/src/SU2_DOT.cpp
index 58d90a193f96..ff4e98c3e297 100644
--- a/SU2_DOT/src/SU2_DOT.cpp
+++ b/SU2_DOT/src/SU2_DOT.cpp
@@ -6,7 +6,7 @@
  *
  * SU2 Project Website: https://su2code.github.io
  *
- * The SU2 Project is maintained by the SU2 Foundation 
+ * The SU2 Project is maintained by the SU2 Foundation
  * (http://su2foundation.org)
  *
  * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md)
@@ -31,16 +31,10 @@ using namespace std;
 
 int main(int argc, char *argv[]) {
 
-  unsigned short iZone, nZone = SINGLE_ZONE, iInst;
+  unsigned short iZone, iInst;
   su2double StartTime = 0.0, StopTime = 0.0, UsedTime = 0.0;
 
   char config_file_name[MAX_STRING_SIZE];
-  ofstream Gradient_file;
-  bool fem_solver = false;
-
-  su2double** Gradient;
-  unsigned short iDV, iDV_Value;
-  int rank, size;
 
   /*--- MPI initialization, and buffer setting ---*/
 
@@ -56,17 +50,17 @@ int main(int argc, char *argv[]) {
   SU2_Comm MPICommunicator(0);
 #endif
 
-  rank = SU2_MPI::GetRank();
-  size = SU2_MPI::GetSize();
+  const int rank = SU2_MPI::GetRank();
+  const int size = SU2_MPI::GetSize();
 
   /*--- Pointer to different structures that will be used throughout the entire code ---*/
 
-  CConfig **config_container            = NULL;
-  CConfig *driver_config                = NULL;
-  CGeometry ***geometry_container       = NULL;
-  CSurfaceMovement **surface_movement   = NULL;
-  CVolumetricMovement **grid_movement   = NULL;
-  unsigned short *nInst                 = NULL;
+  CConfig **config_container            = nullptr;
+  CConfig *driver_config                = nullptr;
+  CGeometry ***geometry_container       = nullptr;
+  CSurfaceMovement **surface_movement   = nullptr;
+  CVolumetricMovement **grid_movement   = nullptr;
+  unsigned short *nInst                 = nullptr;
 
   /*--- Load in the number of zones and spatial dimensions in the mesh file (if no config
    file is specified, default.cfg is used) ---*/
@@ -78,27 +72,23 @@ int main(int argc, char *argv[]) {
    file the number of zones and dimensions from the numerical grid (required
    for variables allocation)  ---*/
 
-  CConfig *config = NULL;
+  CConfig *config = nullptr;
   config = new CConfig(config_file_name, SU2_DOT);
 
-  nZone    = config->GetnZone();
+  const auto nZone = config->GetnZone();
 
   /*--- Definition of the containers per zones ---*/
 
-  config_container    = new CConfig*[nZone];
-  geometry_container  = new CGeometry**[nZone];
-  surface_movement    = new CSurfaceMovement*[nZone];
-  grid_movement       = new CVolumetricMovement*[nZone];
+  config_container    = new CConfig*[nZone] ();
+  geometry_container  = new CGeometry**[nZone] ();
+  surface_movement    = new CSurfaceMovement*[nZone] ();
+  grid_movement       = new CVolumetricMovement*[nZone] ();
 
   nInst               = new unsigned short[nZone];
-  driver_config       = NULL;
+  driver_config       = nullptr;
 
   for (iZone = 0; iZone < nZone; iZone++) {
-    config_container[iZone]       = NULL;
-    geometry_container[iZone]     = NULL;
-    grid_movement [iZone]         = NULL;
-    surface_movement[iZone]       = NULL;
-    nInst[iZone]                  = 1;
+    nInst[iZone] = 1;
   }
 
   /*--- Initialize the configuration of the driver ---*/
@@ -144,13 +134,7 @@ int main(int argc, char *argv[]) {
 
     /*--- Determine whether or not the FEM solver is used, which decides the
      type of geometry classes that are instantiated. ---*/
-    fem_solver = ((config_container[iZone]->GetKind_Solver() == FEM_EULER)          ||
-                  (config_container[iZone]->GetKind_Solver() == FEM_NAVIER_STOKES)  ||
-                  (config_container[iZone]->GetKind_Solver() == FEM_RANS)           ||
-                  (config_container[iZone]->GetKind_Solver() == FEM_LES)            ||
-                  (config_container[iZone]->GetKind_Solver() == DISC_ADJ_FEM_EULER) ||
-                  (config_container[iZone]->GetKind_Solver() == DISC_ADJ_FEM_NS)    ||
-                  (config_container[iZone]->GetKind_Solver() == DISC_ADJ_FEM_RANS));
+    const bool fem_solver = config_container[iZone]->GetFEMSolver();
 
     /*--- Read the number of instances for each zone ---*/
 
@@ -162,7 +146,7 @@ int main(int argc, char *argv[]) {
 
       /*--- Definition of the geometry class to store the primal grid in the partitioning process. ---*/
 
-      CGeometry *geometry_aux = NULL;
+      CGeometry *geometry_aux = nullptr;
 
       /*--- All ranks process the grid and call ParMETIS for partitioning ---*/
 
@@ -318,13 +302,10 @@ int main(int argc, char *argv[]) {
 
       /*--- Initialize structure to store the gradient ---*/
 
-      Gradient = new su2double*[config_container[ZONE_0]->GetnDV()];
+      su2double** Gradient = new su2double*[config_container[ZONE_0]->GetnDV()];
 
-      for (iDV = 0; iDV  < config_container[iZone]->GetnDV(); iDV++){
-        Gradient[iDV] = new su2double[config_container[iZone]->GetnDV_Value(iDV)];
-        for (iDV_Value = 0; iDV_Value < config_container[iZone]->GetnDV_Value(iDV); iDV_Value++){
-          Gradient[iDV][iDV_Value] = 0.0;
-        }
+      for (auto iDV = 0u; iDV  < config_container[iZone]->GetnDV(); iDV++) {
+        Gradient[iDV] = new su2double[config_container[iZone]->GetnDV_Value(iDV)] ();
       }
 
       if (rank == MASTER_NODE)
@@ -332,6 +313,7 @@ int main(int argc, char *argv[]) {
 
       /*--- Write the gradient in a external file ---*/
 
+      ofstream Gradient_file;
       if (rank == MASTER_NODE)
         Gradient_file.open(config_container[iZone]->GetObjFunc_Grad_FileName().c_str(), ios::out);
 
@@ -355,27 +337,24 @@ int main(int argc, char *argv[]) {
 
       OutputGradient(Gradient, config_container[iZone], Gradient_file);
 
-      if (rank == MASTER_NODE)
-        Gradient_file.close();
-
-      for (iDV = 0; iDV  < config_container[iZone]->GetnDV(); iDV++){
+      for (auto iDV = 0u; iDV  < config_container[iZone]->GetnDV(); iDV++){
         delete [] Gradient[iDV];
       }
+      delete [] Gradient;
     }
   }
-  delete [] Gradient;
 
   delete config;
-  config = NULL;
+  config = nullptr;
 
   if (rank == MASTER_NODE)
     cout << "\n------------------------- Solver Postprocessing -------------------------" << endl;
 
-  if (geometry_container != NULL) {
+  if (geometry_container != nullptr) {
     for (iZone = 0; iZone < nZone; iZone++) {
-      if (geometry_container[iZone] != NULL) {
+      if (geometry_container[iZone] != nullptr) {
         for (iInst = 0; iInst < nInst[iZone]; iInst++){
-          if (geometry_container[iZone][iInst] != NULL) {
+          if (geometry_container[iZone][iInst] != nullptr) {
             delete geometry_container[iZone][iInst];
           }
         }
@@ -386,9 +365,9 @@ int main(int argc, char *argv[]) {
   }
   if (rank == MASTER_NODE) cout << "Deleted CGeometry container." << endl;
 
-  if (surface_movement != NULL) {
+  if (surface_movement != nullptr) {
     for (iZone = 0; iZone < nZone; iZone++) {
-      if (surface_movement[iZone] != NULL) {
+      if (surface_movement[iZone] != nullptr) {
         delete surface_movement[iZone];
       }
     }
@@ -396,9 +375,9 @@ int main(int argc, char *argv[]) {
   }
   if (rank == MASTER_NODE) cout << "Deleted CSurfaceMovement class." << endl;
 
-  if (grid_movement != NULL) {
+  if (grid_movement != nullptr) {
     for (iZone = 0; iZone < nZone; iZone++) {
-      if (grid_movement[iZone] != NULL) {
+      if (grid_movement[iZone] != nullptr) {
         delete grid_movement[iZone];
       }
     }
@@ -407,10 +386,10 @@ int main(int argc, char *argv[]) {
   if (rank == MASTER_NODE) cout << "Deleted CVolumetricMovement class." << endl;
 
   delete config;
-  config = NULL;
-  if (config_container != NULL) {
+  config = nullptr;
+  if (config_container != nullptr) {
     for (iZone = 0; iZone < nZone; iZone++) {
-      if (config_container[iZone] != NULL) {
+      if (config_container[iZone] != nullptr) {
         delete config_container[iZone];
       }
     }
@@ -471,7 +450,7 @@ void SetProjection_FD(CGeometry *geometry, CConfig *config, CSurfaceMovement *su
 
   unsigned short nFFDBox = MAX_NUMBER_FFD;
   FFDBox = new CFreeFormDefBox*[nFFDBox];
-  for (iFFDBox = 0; iFFDBox < MAX_NUMBER_FFD; iFFDBox++) FFDBox[iFFDBox] = NULL;
+  for (iFFDBox = 0; iFFDBox < MAX_NUMBER_FFD; iFFDBox++) FFDBox[iFFDBox] = nullptr;
 
   for (iDV = 0; iDV  < nDV; iDV++){
     nDV_Value = config->GetnDV_Value(iDV);
@@ -692,9 +671,9 @@ void SetProjection_FD(CGeometry *geometry, CConfig *config, CSurfaceMovement *su
 
   /*--- Delete memory for parameterization. ---*/
 
-  if (FFDBox != NULL) {
+  if (FFDBox != nullptr) {
     for (iFFDBox = 0; iFFDBox < MAX_NUMBER_FFD; iFFDBox++) {
-      if (FFDBox[iFFDBox] != NULL) {
+      if (FFDBox[iFFDBox] != nullptr) {
         delete FFDBox[iFFDBox];
       }
     }
@@ -718,7 +697,7 @@ void SetProjection_AD(CGeometry *geometry, CConfig *config, CSurfaceMovement *su
   nDim    = geometry->GetnDim();
   nDV     = config->GetnDV();
 
-  VarCoord = NULL;
+  VarCoord = nullptr;
 
   /*--- Discrete adjoint gradient computation ---*/
 
@@ -888,9 +867,8 @@ void SetSensitivity_Files(CGeometry ***geometry, CConfig **config, unsigned shor
 
   unsigned short iZone;
 
-  CSolver *solver  = NULL;
-  COutput *output  = NULL;
-
+  CSolver *solver  = nullptr;
+  COutput *output  = nullptr;
 
   for (iZone = 0; iZone < val_nZone; iZone++) {
 
@@ -1000,4 +978,3 @@ void SetSensitivity_Files(CGeometry ***geometry, CConfig **config, unsigned shor
   }
 
 }
-
diff --git a/SU2_SOL/src/SU2_SOL.cpp b/SU2_SOL/src/SU2_SOL.cpp
index 85b850527aeb..81734122abd5 100644
--- a/SU2_SOL/src/SU2_SOL.cpp
+++ b/SU2_SOL/src/SU2_SOL.cpp
@@ -6,7 +6,7 @@
  *
  * SU2 Project Website: https://su2code.github.io
  *
- * The SU2 Project is maintained by the SU2 Foundation 
+ * The SU2 Project is maintained by the SU2 Foundation
  * (http://su2foundation.org)
  *
  * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md)
@@ -32,14 +32,10 @@ using namespace std;
 
 int main(int argc, char *argv[]) {
 
-  unsigned short iZone, nZone = SINGLE_ZONE, iInst;
+  unsigned short iZone, iInst;
   su2double StartTime = 0.0, StopTime = 0.0, UsedTime = 0.0;
 
   char config_file_name[MAX_STRING_SIZE];
-  int rank = MASTER_NODE;
-  int size = SINGLE_NODE;
-  bool fem_solver = false;
-  bool multizone = false;
 
   /*--- MPI initialization ---*/
 
@@ -50,12 +46,12 @@ int main(int argc, char *argv[]) {
   SU2_Comm MPICommunicator(0);
 #endif
 
-  rank = SU2_MPI::GetRank();
-  size = SU2_MPI::GetSize();
+  const int rank = SU2_MPI::GetRank();
+  const int size = SU2_MPI::GetSize();
 
   /*--- Pointer to different structures that will be used throughout the entire code ---*/
 
-  COutput **output                 = NULL;
+  COutput **output                = NULL;
   CGeometry ***geometry_container = NULL;
   CSolver ***solver_container     = NULL;
   CConfig **config_container      = NULL;
@@ -71,23 +67,19 @@ int main(int argc, char *argv[]) {
   CConfig *config = NULL;
   config = new CConfig(config_file_name, SU2_SOL);
 
-  nZone = config->GetnZone();
+  const auto nZone = config->GetnZone();
 
   /*--- Definition of the containers per zones ---*/
 
-  solver_container = new CSolver**[nZone];
-  config_container = new CConfig*[nZone];
-  geometry_container = new CGeometry**[nZone];
+  solver_container = new CSolver**[nZone] ();
+  config_container = new CConfig*[nZone] ();
+  geometry_container = new CGeometry**[nZone] ();
   nInst = new unsigned short[nZone];
   driver_config = NULL;
-  output   = new COutput*[nZone];
-  
+  output = new COutput*[nZone] ();
+
   for (iZone = 0; iZone < nZone; iZone++) {
-    solver_container[iZone]       = NULL;
-    config_container[iZone]       = NULL;
-    geometry_container[iZone]     = NULL;
-    nInst[iZone]                  = 1;
-    output[iZone]                 = NULL;
+    nInst[iZone] = 1;
   }
 
   /*--- Initialize the configuration of the driver ---*/
@@ -97,7 +89,7 @@ int main(int argc, char *argv[]) {
   char zone_file_name[MAX_STRING_SIZE];
 
   /*--- Store a boolean for multizone problems ---*/
-  multizone = (config->GetMultizone_Problem());
+  const bool multizone = config->GetMultizone_Problem();
 
   /*--- Loop over all zones to initialize the various classes. In most
    cases, nZone is equal to one. This represents the solution of a partial
@@ -122,7 +114,7 @@ int main(int argc, char *argv[]) {
   }
 
   /*--- Set the multizone part of the problem. ---*/
-  if (config->GetMultizone_Problem()){
+  if (multizone) {
     for (iZone = 0; iZone < nZone; iZone++) {
       /*--- Set the interface markers for multizone ---*/
       config_container[iZone]->SetMultizone(driver_config, config_container);
@@ -134,13 +126,7 @@ int main(int argc, char *argv[]) {
 
     /*--- Determine whether or not the FEM solver is used, which decides the
      type of geometry classes that are instantiated. ---*/
-    fem_solver = ((config_container[iZone]->GetKind_Solver() == FEM_EULER)          ||
-                  (config_container[iZone]->GetKind_Solver() == FEM_NAVIER_STOKES)  ||
-                  (config_container[iZone]->GetKind_Solver() == FEM_RANS)           ||
-                  (config_container[iZone]->GetKind_Solver() == FEM_LES)            ||
-                  (config_container[iZone]->GetKind_Solver() == DISC_ADJ_FEM_EULER) ||
-                  (config_container[iZone]->GetKind_Solver() == DISC_ADJ_FEM_NS)    ||
-                  (config_container[iZone]->GetKind_Solver() == DISC_ADJ_FEM_RANS));
+    const bool fem_solver = config_container[iZone]->GetFEMSolver();
 
     /*--- Read the number of instances for each zone ---*/
 
@@ -169,8 +155,7 @@ int main(int argc, char *argv[]) {
       if ( fem_solver ) geometry_aux->SetColorFEMGrid_Parallel(config_container[iZone]);
       else              geometry_aux->SetColorGrid_Parallel(config_container[iZone]);
 
-      /*--- Allocate the memory of the current domain, and
-     divide the grid between the nodes ---*/
+      /*--- Allocate the memory of the current domain, and divide the grid between the nodes ---*/
 
       geometry_container[iZone][iInst] = NULL;
 
@@ -195,7 +180,7 @@ int main(int argc, char *argv[]) {
       /*--- Add the Send/Receive boundaries ---*/
 
       geometry_container[iZone][iInst]->SetSendReceive(config_container[iZone]);
-      
+
       /*--- Add the Send/Receive boundaries ---*/
 
       geometry_container[iZone][iInst]->SetBoundaries(config_container[iZone]);
@@ -211,9 +196,9 @@ int main(int argc, char *argv[]) {
       geometry_container[iZone][iInst]->SetGlobal_to_Local_Point();
 
       /*--- Create the point-to-point MPI communication structures for the fvm solver. ---*/
-      
+
       if (!fem_solver) geometry_container[iZone][iInst]->PreprocessP2PComms(geometry_container[iZone][iInst], config_container[iZone]);
-      
+
       /* Test for a fem solver, because some more work must be done. */
 
       if (fem_solver) {
@@ -235,9 +220,8 @@ int main(int argc, char *argv[]) {
 
   }
 
-  /*--- Determine whether the simulation is a FSI simulation ---*/
-
-  bool fsi = config_container[ZONE_0]->GetFSI_Simulation();
+  const bool fsi = config_container[ZONE_0]->GetFSI_Simulation();
+  const bool fem_solver = config_container[ZONE_0]->GetFEMSolver();
 
   /*--- Set up a timer for performance benchmarking (preprocessing time is included) ---*/
 
@@ -249,13 +233,11 @@ int main(int argc, char *argv[]) {
 
   if (rank == MASTER_NODE)
     cout << endl <<"------------------------- Solution Postprocessing -----------------------" << endl;
-  
+
   /*---  Check whether this is an FSI, fluid unsteady, harmonic balance or structural dynamic simulation and call the
    solution merging routines accordingly.---*/
 
-  if (multizone){
-    
-
+  if (multizone) {
 
     bool TimeDomain = driver_config->GetTime_Domain();
 
@@ -278,11 +260,11 @@ int main(int argc, char *argv[]) {
         config_container[iZone]->SetiInst(INST_0);
         config_container[iZone]->SetTimeIter(TimeIter);
         solver_container[iZone][INST_0] = new CBaselineSolver(geometry_container[iZone][INST_0], config_container[iZone]);
-        
+
         output[iZone] = new CBaselineOutput(config_container[iZone], geometry_container[iZone][INST_0]->GetnDim(), solver_container[iZone][INST_0]);
         output[iZone]->PreprocessVolumeOutput(config_container[iZone]);
         output[iZone]->PreprocessHistoryOutput(config_container[iZone], false);
-        
+
       }
 
       /*--- Loop over the whole time domain ---*/
@@ -309,11 +291,11 @@ int main(int argc, char *argv[]) {
             solver_container[iZone][INST_0]->LoadRestart(geometry_container[iZone], &solver_container[iZone], config_container[iZone], TimeIter, true);
           }
 
-          
+
           for (iZone = 0; iZone < nZone; iZone++){
-            
+
             WriteFiles(config_container[iZone], geometry_container[iZone][INST_0], &solver_container[iZone][INST_0], output[iZone], TimeIter);
-            
+
           }
         }
 
@@ -332,12 +314,12 @@ int main(int argc, char *argv[]) {
         output[iZone] = new CBaselineOutput(config_container[iZone], geometry_container[iZone][INST_0]->GetnDim(), solver_container[iZone][INST_0]);
         output[iZone]->PreprocessVolumeOutput(config_container[iZone]);
         output[iZone]->PreprocessHistoryOutput(config_container[iZone], false);
-      
+
       }
       for (iZone = 0; iZone < nZone; iZone++){
-        
-        WriteFiles(config_container[iZone], geometry_container[iZone][INST_0], &solver_container[iZone][INST_0], output[iZone], 0);      
-        
+
+        WriteFiles(config_container[iZone], geometry_container[iZone][INST_0], &solver_container[iZone][INST_0], output[iZone], 0);
+
       }
     }
 
@@ -410,7 +392,7 @@ int main(int argc, char *argv[]) {
           output[ZONE_0] = new CBaselineOutput(config_container[ZONE_0], geometry_container[ZONE_0][INST_0]->GetnDim(), solver_container[ZONE_0][INST_0]);
           output[ZONE_0]->PreprocessVolumeOutput(config_container[ZONE_0]);
           output[ZONE_0]->PreprocessHistoryOutput(config_container[ZONE_0], false);
-         
+
           SolutionInstantiatedFlow = true;
         }
           solver_container[ZONE_0][INST_0]->LoadRestart_FSI(geometry_container[ZONE_0][INST_0], config_container[ZONE_0], TimeIter);
@@ -432,9 +414,9 @@ int main(int argc, char *argv[]) {
 
         if (rank == MASTER_NODE) cout << "Writing the volume solution for time step " << TimeIter << "." << endl;
         for (iZone = 0; iZone < nZone; iZone++){
-          
+
           WriteFiles(config_container[iZone], geometry_container[iZone][INST_0], &solver_container[iZone][INST_0], output[iZone], TimeIter);
-          
+
         }
       }
 
@@ -501,18 +483,18 @@ int main(int argc, char *argv[]) {
 
               if (rank == MASTER_NODE)
                 cout << "Writing the volume solution for time step " << TimeIter << "." << endl;
-              
+
               for (iZone = 0; iZone < nZone; iZone++){
-                
+
                 WriteFiles(config_container[iZone], geometry_container[iZone][INST_0], &solver_container[iZone][INST_0], output[iZone], TimeIter);
-                
+
               }
         }
-        
+
         TimeIter++;
         if (StopCalc) break;
       }
-      
+
     } else {
 
     /*--- Steady simulation: merge the single solution file. ---*/
@@ -526,14 +508,14 @@ int main(int argc, char *argv[]) {
       output[iZone]->PreprocessHistoryOutput(config_container[ZONE_0], false);
       solver_container[iZone][INST_0]->LoadRestart(&geometry_container[iZone][INST_0], &solver_container[iZone], config_container[iZone], 0, true);
     }
-    
+
     for (iZone = 0; iZone < nZone; iZone++){
-      
+
       WriteFiles(config_container[iZone], geometry_container[iZone][INST_0], &solver_container[iZone][INST_0], output[iZone], 0);
-      
-    }   
+
+    }
     }
-    
+
   }
   else {
 
@@ -588,23 +570,23 @@ int main(int argc, char *argv[]) {
                   output[iZone] = new CBaselineOutput(config_container[iZone], geometry_container[iZone][INST_0]->GetnDim(), solver_container[iZone][INST_0]);
                   output[iZone]->PreprocessVolumeOutput(config_container[iZone]);
                   output[iZone]->PreprocessHistoryOutput(config_container[iZone], false);
-                  
+
                   SolutionInstantiated[iZone] = true;
                 }
                   config_container[iZone]->SetiInst(INST_0);
                   solver_container[iZone][INST_0]->LoadRestart(geometry_container[iZone], &solver_container[iZone], config_container[iZone], TimeIter, true);
               }
-              
+
               if (rank == MASTER_NODE)
                 cout << "Writing the volume solution for time step " << TimeIter << "." << endl;
-              
+
               for (iZone = 0; iZone < nZone; iZone++){
-                
+
                 WriteFiles(config_container[iZone], geometry_container[iZone][INST_0], &solver_container[iZone][INST_0], output[iZone], TimeIter);
-                
-              } 
-              
-              
+
+              }
+
+
         }
 
         TimeIter++;
@@ -626,7 +608,7 @@ int main(int argc, char *argv[]) {
           /*--- Either instantiate the solution class or load a restart file. ---*/
           solver_container[iZone][iInst] = new CBaselineSolver(geometry_container[iZone][iInst], config_container[iZone]);
           solver_container[iZone][iInst]->LoadRestart(geometry_container[iZone], &solver_container[iZone], config_container[iZone], iInst, true);
-          output[iZone] = new CBaselineOutput(config_container[iZone], geometry_container[iZone][iInst]->GetnDim(), solver_container[iZone][iInst]);          
+          output[iZone] = new CBaselineOutput(config_container[iZone], geometry_container[iZone][iInst]->GetnDim(), solver_container[iZone][iInst]);
           output[iZone]->PreprocessVolumeOutput(config_container[iZone]);
           output[iZone]->PreprocessHistoryOutput(config_container[iZone], false);
 
@@ -636,7 +618,7 @@ int main(int argc, char *argv[]) {
           }
 
           WriteFiles(config_container[iZone], geometry_container[iZone][iInst], &solver_container[iZone][iInst], output[iZone], iInst);
-          
+
         }
 
       }
@@ -691,7 +673,7 @@ int main(int argc, char *argv[]) {
                   output[iZone] = new CBaselineOutput(config_container[iZone], geometry_container[iZone][INST_0]->GetnDim(), solver_container[iZone][INST_0]);
                   output[iZone]->PreprocessVolumeOutput(config_container[iZone]);
                   output[iZone]->PreprocessHistoryOutput(config_container[iZone], false);
-                  
+
                   SolutionInstantiated = true;
                 }
                 config_container[iZone]->SetiInst(INST_0);
@@ -701,22 +683,22 @@ int main(int argc, char *argv[]) {
               if (rank == MASTER_NODE)
                 cout << "Writing the volume solution for time step " << TimeIter << "." << endl;
               for (iZone = 0; iZone < nZone; iZone++){
-                
+
                 WriteFiles(config_container[iZone], geometry_container[iZone][INST_0], &solver_container[iZone][INST_0], output[iZone], TimeIter);
-                
-              }        
+
+              }
         }
-        
+
         TimeIter++;
         if (StopCalc) break;
       }
-      
+
 		  }
-    
+
     else {
 
       /*--- Steady simulation: merge the single solution file. ---*/
-      
+
       for (iZone = 0; iZone < nZone; iZone++) {
         config_container[iZone]->SetiInst(INST_0);
         /*--- Definition of the solution class ---*/
@@ -725,23 +707,23 @@ int main(int argc, char *argv[]) {
         output[iZone] = new CBaselineOutput(config_container[iZone], geometry_container[iZone][INST_0]->GetnDim(), solver_container[iZone][INST_0]);
         output[iZone]->PreprocessVolumeOutput(config_container[iZone]);
         output[iZone]->PreprocessHistoryOutput(config_container[iZone], false);
-        
+
       }
       for (iZone = 0; iZone < nZone; iZone++){
-        
+
         WriteFiles(config_container[iZone], geometry_container[iZone][INST_0], &solver_container[iZone][INST_0], output[iZone], 0);
-        
-      }  
+
+      }
     }
-    
+
   }
-  
+
   delete config;
   config = NULL;
 
   if (rank == MASTER_NODE)
     cout << endl <<"------------------------- Solver Postprocessing -------------------------" << endl;
-  
+
   if (geometry_container != NULL) {
     for (iZone = 0; iZone < nZone; iZone++) {
       for (iInst = 0; iInst < nInst[iZone]; iInst++){
@@ -755,7 +737,7 @@ int main(int argc, char *argv[]) {
     delete [] geometry_container;
   }
   if (rank == MASTER_NODE) cout << "Deleted CGeometry container." << endl;
-  
+
   if (solver_container != NULL) {
     for (iZone = 0; iZone < nZone; iZone++) {
       for (iInst = 0; iInst < nInst[iZone]; iInst++){
@@ -769,7 +751,7 @@ int main(int argc, char *argv[]) {
     delete [] solver_container;
   }
   if (rank == MASTER_NODE) cout << "Deleted CSolver class." << endl;
-  
+
   if (config_container != NULL) {
     for (iZone = 0; iZone < nZone; iZone++) {
       if (config_container[iZone] != NULL) {
@@ -779,7 +761,7 @@ int main(int argc, char *argv[]) {
     delete [] config_container;
   }
   if (rank == MASTER_NODE) cout << "Deleted CConfig container." << endl;
-  
+
   if (output != NULL) {
     for (iZone = 0; iZone < nZone; iZone++) {
       if (output[iZone] != NULL) {
@@ -789,54 +771,54 @@ int main(int argc, char *argv[]) {
     delete [] output;
   }
   if (rank == MASTER_NODE) cout << "Deleted COutput class." << endl;
-  
+
   /*--- Synchronization point after a single solver iteration. Compute the
    wall clock time required. ---*/
-  
+
 #ifdef HAVE_MPI
   StopTime = MPI_Wtime();
 #else
   StopTime = su2double(clock())/su2double(CLOCKS_PER_SEC);
 #endif
-  
+
   /*--- Compute/print the total time for performance benchmarking. ---*/
-  
+
   UsedTime = StopTime-StartTime;
   if (rank == MASTER_NODE) {
     cout << "\nCompleted in " << fixed << UsedTime << " seconds on "<< size;
     if (size == 1) cout << " core." << endl; else cout << " cores." << endl;
   }
-  
+
   /*--- Exit the solver cleanly ---*/
-  
+
   if (rank == MASTER_NODE)
     cout << endl <<"------------------------- Exit Success (SU2_SOL) ------------------------" << endl << endl;
-  
+
   /*--- Finalize MPI parallelization ---*/
-  
+
 #ifdef HAVE_MPI
   SU2_MPI::Finalize();
 #endif
-  
+
   return EXIT_SUCCESS;
 }
 
 void WriteFiles(CConfig *config, CGeometry* geometry, CSolver** solver_container, COutput *output, unsigned long TimeIter){
-  
+
   /*--- Load history data (volume output might require some values) --- */
-  
+
   output->SetHistory_Output(geometry, solver_container, config, TimeIter, 0, 0);
-  
+
   /*--- Load the data --- */
-  
+
   output->Load_Data(geometry, config, solver_container);
-  
+
   /*--- Set the filenames ---*/
-  
+
   output->SetVolume_Filename(config->GetVolume_FileName());
-  
+
   output->SetSurface_Filename(config->GetSurfCoeff_FileName());
-  
+
   for (unsigned short iFile = 0; iFile < config->GetnVolumeOutputFiles(); iFile++){
     unsigned short* FileFormat = config->GetVolumeOutputFiles();
     if (FileFormat[iFile] != RESTART_ASCII &&
@@ -844,5 +826,5 @@ void WriteFiles(CConfig *config, CGeometry* geometry, CSolver** solver_container
         FileFormat[iFile] != CSV)
       output->WriteToFile(config, geometry, FileFormat[iFile]);
   }
-  
+
 }