diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 794622bb1a1e..6f837fd0bc54 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -1277,6 +1277,18 @@ class CConfig {
   /*--- Additional flamelet solver options ---*/
   FluidFlamelet_ParsedOptions flamelet_ParsedOptions; /*!< \brief Additional flamelet solver options */
 
+  /*--- Mesh adaptation options ---*/
+  bool Compute_Metric;                     /*!< \brief Determines if error estimation is taking place */
+  bool Normalize_Metric;                   /*!< \brief Determines if metric tensor normalization is taking place */
+  unsigned short Kind_Hessian_Method;      /*!< \brief Numerical method for computation of Hessians. */
+  unsigned short nMetric_Sensor;           /*!< \brief Number of sensors to use for adaptation. */
+  METRIC_SENSOR* Metric_Sensor;            /*!< \brief Sensors to use for adaptation. */
+  unsigned short Metric_Norm;              /*!< \brief Lp-norm for mesh adaptation */
+  unsigned long Metric_Complexity;         /*!< \brief Constraint mesh complexity */
+  su2double Metric_Hmax,                   /*!< \brief Maximum cell size */
+            Metric_Hmin,                   /*!< \brief Minimum cell size */
+            Metric_ARmax;                  /*!< \brief Maximum cell aspect ratio */
+
   /*!
    * \brief Set the default values of config options not set in the config file using another config object.
    * \param config - Config object to use the default values from.
@@ -10142,4 +10154,101 @@ class CConfig {
    */
   const FluidFlamelet_ParsedOptions& GetFlameletParsedOptions() const { return flamelet_ParsedOptions; }
 
+  /*!
+   * \brief Check if a metric tensor field is being computed
+   * \return <code>TRUE<\code> if a metric tensor field is being computed
+   */
+  bool GetCompute_Metric(void) const { return Compute_Metric; }
+
+  /*!
+   * \brief Check if metric tensor normalization is being carried out
+   * \return <code>TRUE<\code> if metric normalization is taking place
+   */
+  bool GetNormalize_Metric(void) const { return Normalize_Metric; }
+
+  /*!
+   * \brief Get the kind of method for computation of Hessians used for anisotropy.
+   * \return Numerical method for computation of Hessians used for anisotropy.
+   */
+  unsigned short GetKind_Hessian_Method(void) const { return Kind_Hessian_Method; }
+
+  /*!
+   * \brief Get adaptation sensor
+   */
+  METRIC_SENSOR GetMetric_Sensor(unsigned short iSens) const { return Metric_Sensor[iSens]; }
+
+  /*!
+   * \brief Get corresponding string from metric sensor type
+   */
+  string GetMetric_SensorString(unsigned short iSens) const {
+    string sensor_name;
+    switch (Metric_Sensor[iSens]) {
+      case METRIC_SENSOR::DENSITY:
+        sensor_name = "Density";
+        break;
+      case METRIC_SENSOR::MACH:
+        sensor_name = "Mach";
+        break;
+      case METRIC_SENSOR::PRESSURE:
+        sensor_name = "Pressure";
+        break;
+      case METRIC_SENSOR::TOTAL_PRESSURE:
+        sensor_name = "Total Pressure";
+        break;
+      case METRIC_SENSOR::TEMPERATURE:
+        sensor_name = "Temperature";
+        break;
+      case METRIC_SENSOR::TEMPERATURE_VE:
+        sensor_name = "Temperature_ve";
+        break;
+      case METRIC_SENSOR::ENERGY:
+        sensor_name = "Energy";
+        break;
+      case METRIC_SENSOR::ENERGY_VE:
+        sensor_name = "Energy_ve";
+        break;
+      case METRIC_SENSOR::GOAL:
+        sensor_name = "Goal";
+        break;
+      default:
+        SU2_MPI::Error("Unsupported metric sensor.", CURRENT_FUNCTION);
+    }
+
+    return sensor_name;
+  }
+
+  /*!
+   * \brief Get number of adaptation sensors
+   */
+  unsigned short GetnMetric_Sensor(void) const { return nMetric_Sensor; }
+
+  /*!
+   * \brief Get adaptation norm value (Lp)
+   */
+  unsigned short GetMetric_Norm(void) const { return Metric_Norm; }
+
+  /*!
+   * \brief Get maximum cell size
+   * \return Maximum cell size
+   */
+  su2double GetMetric_Hmax(void) const { return Metric_Hmax; }
+
+  /*!
+   * \brief Get minimum cell size
+   * \return Minimum cell size
+   */
+  su2double GetMetric_Hmin(void) const { return Metric_Hmin; }
+
+  /*!
+   * \brief Get maximum cell aspect ratio
+   * \return Maximum cell aspect ratio
+   */
+  su2double GetMetric_ARmax(void) const { return Metric_ARmax; }
+
+  /*!
+   * \brief Get constraint complexity
+   * \return Mesh complexity
+   */
+  unsigned long GetMetric_Complexity(void) const { return Metric_Complexity; }
+
 };
diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index bf99257b1387..6b7320b62bed 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -2621,6 +2621,8 @@ enum PERIODIC_QUANTITIES {
   PERIODIC_LIM_PRIM_1 ,  /*!< \brief Primitive limiter communication phase 1 of 2 (periodic only). */
   PERIODIC_LIM_PRIM_2 ,  /*!< \brief Primitive limiter communication phase 2 of 2 (periodic only). */
   PERIODIC_IMPLICIT   ,  /*!< \brief Implicit update communication to ensure consistency across periodic boundaries. */
+  PERIODIC_GRAD_ADAPT ,  /*!< \brief Gradient vectors for anisotropic sizing metric (periodic only). */
+  PERIODIC_HESSIAN    ,  /*!< \brief Hessian tensors for anisotropic sizing metric (periodic only). */
 };
 
 /*!
@@ -2652,6 +2654,8 @@ enum class MPI_QUANTITIES {
   MESH_DISPLACEMENTS   ,  /*!< \brief Mesh displacements at the interface. */
   SOLUTION_TIME_N      ,  /*!< \brief Solution at time n. */
   SOLUTION_TIME_N1     ,  /*!< \brief Solution at time n-1. */
+  GRADIENT_ADAPT       ,  /*!< \brief Gradient vectors for anisotropic sizing metric tensor. */
+  HESSIAN              ,  /*!< \brief Hessian tensors for anisotropic sizing metric tensor. */
 };
 
 /*!
@@ -2797,6 +2801,32 @@ static const MapType<std::string, ENUM_SOBOLEV_MODUS> Sobolev_Modus_Map = {
   MakePair("ONLY_GRADIENT",        ENUM_SOBOLEV_MODUS::ONLY_GRAD)
 };
 
+/*!
+ * \brief Types of sensors for anisotropic metric
+ */
+enum class METRIC_SENSOR {
+  DENSITY,         /*!< \brief Density feature-based metric. */
+  MACH,            /*!< \brief Mach feature-based metric. */
+  PRESSURE,        /*!< \brief Pressure feature-based metric. */
+  TOTAL_PRESSURE,  /*!< \brief Total pressure feature-based metric. */
+  TEMPERATURE,     /*!< \brief Temperature feature-based metric. */
+  TEMPERATURE_VE,  /*!< \brief Vibrational/electronic temperature feature-based metric. */
+  ENERGY,          /*!< \brief Energy feature-based metric. */
+  ENERGY_VE,       /*!< \brief Vibrational/electronic energy feature-based metric. */
+  GOAL,            /*!< \brief Goal-oriented metric. */
+};
+static const MapType<std::string, METRIC_SENSOR> Metric_Sensor_Map = {
+  MakePair("DENSITY", METRIC_SENSOR::DENSITY)
+  MakePair("MACH", METRIC_SENSOR::MACH)
+  MakePair("PRESSURE", METRIC_SENSOR::PRESSURE)
+  MakePair("TOTAL_PRESSURE", METRIC_SENSOR::TOTAL_PRESSURE)
+  MakePair("TEMPERATURE", METRIC_SENSOR::TEMPERATURE)
+  MakePair("TEMPERATURE_VE", METRIC_SENSOR::TEMPERATURE_VE)
+  MakePair("ENERGY", METRIC_SENSOR::ENERGY)
+  MakePair("ENERGY_VE", METRIC_SENSOR::ENERGY_VE)
+  MakePair("GOAL", METRIC_SENSOR::GOAL)
+};
+
 /*!
  * \brief Type of mesh deformation
  */
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 0a9cca0f63f7..956d6bcf4653 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -3034,6 +3034,55 @@ void CConfig::SetConfig_Options() {
   /*!\brief ROM_SAVE_FREQ \n DESCRIPTION: How often to save snapshots for unsteady problems.*/
   addUnsignedShortOption("ROM_SAVE_FREQ", rom_save_freq, 1);
 
+  /*--- options that are used for mesh adaptation ---*/
+  /*!\par CONFIG_CATEGORY:Adaptation Options \ingroup Config*/
+
+  /*!\brief COMPUTE_METRIC \n DESCRIPTION: Compute a metric tensor field */
+  addBoolOption("COMPUTE_METRIC", Compute_Metric, false);
+  /*!\brief NORMALIZE_METRIC \n DESCRIPTION: Normalize the metric tensor */
+  addBoolOption("NORMALIZE_METRIC", Normalize_Metric, true);
+  /*!\brief NUM_METHOD_HESS \n DESCRIPTION: Numerical method for Hessian computation \n OPTIONS: See \link Gradient_Map \endlink. \n DEFAULT: GREEN_GAUSS. \ingroup Config*/
+  addEnumOption("NUM_METHOD_HESS", Kind_Hessian_Method, Gradient_Map, GREEN_GAUSS);
+
+  /*!\brief METRIC_SENSOR \n DESCRIPTION: Sensors for mesh adaptation */
+  addEnumListOption("METRIC_SENSOR", nMetric_Sensor, Metric_Sensor, Metric_Sensor_Map);
+  /*!\brief METRIC_NORM \n DESCRIPTION: Lp-norm for mesh adaptation */
+  addUnsignedShortOption("METRIC_NORM", Metric_Norm, 2);
+  /*!\brief METRIC_COMPLEXITY \n DESCRIPTION: Constraint mesh complexity */
+  addUnsignedLongOption("METRIC_COMPLEXITY", Metric_Complexity, 10000);
+
+  /*!\brief METRIC_HMAX \n DESCRIPTION: Constraint maximum cell size */
+  addDoubleOption("METRIC_HMAX", Metric_Hmax, 10.0);
+  /*!\brief METRIC_HMIN \n DESCRIPTION: Constraint minimum cell size */
+  addDoubleOption("METRIC_HMIN", Metric_Hmin, 1.0E-8);
+  /*!\brief METRIC_ARMAX \n DESCRIPTION: Constraint maximum cell aspect ratio */
+  addDoubleOption("METRIC_ARMAX", Metric_ARmax, 1.0E6);
+
+  /*!\brief ADAP_ITER \n DESCRIPTION: Mesh adaptation inner iterations per complexity */
+  addPythonOption("ADAP_ITER");
+  /*!\brief ADAP_COMPLEXITIES \n DESCRIPTION: List of constraint (target) mesh complexities for mesh convergence study */
+  addPythonOption("ADAP_COMPLEXITIES");
+  /*!\brief ADAP_FLOW_ITERS \n DESCRIPTION: Primal solver iterations at each target complexity */
+  addPythonOption("ADAP_FLOW_ITERS");
+  /*!\brief ADAP_ADJ_ITERS \n DESCRIPTION: Adjoint solver iterations at each target complexity */
+  addPythonOption("ADAP_ADJ_ITERS");
+  /*!\brief ADAP_FLOW_CFLS \n DESCRIPTION: Primal solver CFL number at each target complexity */
+  addPythonOption("ADAP_FLOW_CFLS");
+  /*!\brief ADAP_ADJ_CFLS \n DESCRIPTION: Adjoint solver CFL number at each target complexity */
+  addPythonOption("ADAP_ADJ_CFLS");
+  /*!\brief ADAP_RESIDUAL_REDUCTIONS \n DESCRIPTION: Residual reduction at each target complexity */
+  addPythonOption("ADAP_RESIDUAL_REDUCTIONS");
+  /*!\brief ADAP_HMAXS \n DESCRIPTION: Maximum cell size at each target complexity */
+  addPythonOption("ADAP_HMAXS");
+  /*!\brief ADAP_HMINS \n DESCRIPTION: Minimum cell size at each target complexity */
+  addPythonOption("ADAP_HMINS");
+  /*!\brief ADAP_HGRAD \n DESCRIPTION: Size gradation smoothing parameter */
+  addPythonOption("ADAP_HGRAD");
+  /*!\brief ADAP_HAUSD \n DESCRIPTION: Hausdorff distance parameter for surface remeshing */
+  addPythonOption("ADAP_HAUSD");
+  /*!\brief ADAP_ANGLE \n DESCRIPTION: Sharp angle detection parameter for surface remeshing */
+  addPythonOption("ADAP_ANGLE");
+
   /* END_CONFIG_OPTIONS */
 
 }
@@ -5752,6 +5801,36 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     SU2_MPI::Error("BOUNDED_SCALAR discretization can only be used for incompressible problems.", CURRENT_FUNCTION);
   }
 
+  /*--- Checks for mesh adaptation ---*/
+  if (Compute_Metric) {
+    /*--- Check that config is valid for requested sensor ---*/
+    for (auto iSensor = 0; iSensor < nMetric_Sensor; iSensor++) {
+      /*--- TODO: If using GOAL, it must be the only sensor and the discrete adjoint must be used ---*/
+      /*--- For now, goal-oriented adaptation is unsupported ---*/
+      if (Metric_Sensor[iSensor] == METRIC_SENSOR::GOAL) {
+        SU2_MPI::Error("Adaptation sensor GOAL not yet supported.", CURRENT_FUNCTION);
+      }
+
+      if (Kind_Solver == MAIN_SOLVER::NEMO_EULER || Kind_Solver == MAIN_SOLVER::NEMO_NAVIER_STOKES) {
+          if (Metric_Sensor[iSensor] == METRIC_SENSOR::DENSITY || Metric_Sensor[iSensor] == METRIC_SENSOR::TOTAL_PRESSURE)
+            SU2_MPI::Error(string("Adaptation sensor ") + GetMetric_SensorString(iSensor) + string(" not available for NEMO problems."), CURRENT_FUNCTION);
+        }
+        if (Kind_Solver != MAIN_SOLVER::NEMO_EULER && Kind_Solver != MAIN_SOLVER::NEMO_NAVIER_STOKES) {
+          if (Metric_Sensor[iSensor] == METRIC_SENSOR::TEMPERATURE_VE || Metric_Sensor[iSensor] == METRIC_SENSOR::ENERGY_VE)
+            SU2_MPI::Error(string("Adaptation sensor ") + GetMetric_SensorString(iSensor) + string(" not available for non-NEMO problems."), CURRENT_FUNCTION);
+        }
+        if (Kind_Solver == MAIN_SOLVER::INC_EULER || Kind_Solver == MAIN_SOLVER::INC_NAVIER_STOKES || Kind_Solver == MAIN_SOLVER::INC_RANS) {
+          if (Metric_Sensor[iSensor] == METRIC_SENSOR::MACH)
+            SU2_MPI::Error(string("Adaptation sensor ") + GetMetric_SensorString(iSensor) + string(" not available for INC problems."), CURRENT_FUNCTION);
+        }
+    }
+
+    /*--- Only GG Hessians for now ---*/
+    if (Kind_Hessian_Method != GREEN_GAUSS) {
+      SU2_MPI::Error("NUM_METHOD_HESS must be GREEN_GAUSS.", CURRENT_FUNCTION);
+    }
+  }
+
 }
 
 void CConfig::SetMarkers(SU2_COMPONENT val_software) {
@@ -7925,6 +8004,30 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
       cout << "Actuator disk BEM method propeller data read from file: " << GetBEM_prop_filename() << endl;
     }
   }
+
+  if (val_software == SU2_COMPONENT::SU2_CFD || val_software == SU2_COMPONENT::SU2_SOL) {
+    if (Compute_Metric) {
+      cout << endl <<"---------------- Mesh Adaptation Information ( Zone "  << iZone << " ) -----------------" << endl;
+      cout << "Adaptation sensor(s): ";
+      for (auto iSensor = 0; iSensor < nMetric_Sensor; iSensor++) {
+        cout << GetMetric_SensorString(iSensor);
+        if (iSensor < nMetric_Sensor - 1 ) cout << ", ";
+      }
+      cout << endl;
+      switch (Kind_Hessian_Method) {
+        case GREEN_GAUSS: cout << "Hessian for adaptive metric: Green-Gauss." << endl; break;
+      }
+      if (Normalize_Metric) {
+        cout << "Target complexity: " << Metric_Complexity << endl;
+        cout << "Lp norm: " << Metric_Norm << endl;
+        cout << "Min. edge length: " << Metric_Hmin << endl;
+        cout << "Max. edge length: " << Metric_Hmax << endl;
+      }
+      else {
+        cout << "Output unnormalized metric field." << endl;
+      }
+    }
+  }
 }
 
 bool CConfig::TokenizeString(string & str, string & option_name,
diff --git a/SU2_CFD/include/drivers/CSinglezoneDriver.hpp b/SU2_CFD/include/drivers/CSinglezoneDriver.hpp
index 355369556ea5..d7e789ae66ce 100644
--- a/SU2_CFD/include/drivers/CSinglezoneDriver.hpp
+++ b/SU2_CFD/include/drivers/CSinglezoneDriver.hpp
@@ -110,4 +110,9 @@ class CSinglezoneDriver : public CDriver {
    */
   bool Monitor(unsigned long TimeIter) override;
 
+  /*!
+   * \brief Perform all steps to compute the metric tensor.
+   */
+  virtual void ComputeMetricField();
+
 };
diff --git a/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp b/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp
index 98a170cf3810..a2f984e1560a 100644
--- a/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp
+++ b/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp
@@ -179,6 +179,123 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER
   solver->InitiateComms(&geometry, &config, kindMpiComm);
   solver->CompleteComms(&geometry, &config, kindMpiComm);
 }
+
+template<size_t nDim, class GradientType>
+void computeHessiansGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PERIODIC_QUANTITIES kindPeriodicComm,
+                               CGeometry& geometry, const CConfig& config, const GradientType& gradient,
+                               const size_t varBegin, const size_t varEnd, const int idxVel, GradientType& hessian) {
+  const size_t nPointDomain = geometry.GetnPointDomain();
+
+#ifdef HAVE_OMP
+  constexpr size_t OMP_MAX_CHUNK = 512;
+
+  const auto chunkSize = computeStaticChunkSize(nPointDomain, omp_get_max_threads(), OMP_MAX_CHUNK);
+#endif
+
+  const size_t nSymMat = 3 * (nDim - 1);
+  su2double diagScale;
+
+  /*--- For each (non-halo) volume integrate over its faces (edges). ---*/
+
+  SU2_OMP_FOR_DYN(chunkSize)
+  for (size_t iPoint = 0; iPoint < nPointDomain; ++iPoint) {
+    auto nodes = geometry.nodes;
+
+    /*--- Clear the Hessian. --*/
+
+    for (size_t iVar = varBegin; iVar < varEnd; ++iVar)
+      for (size_t iMat = 0; iMat < nSymMat; ++iMat) hessian(iPoint, iVar, iMat) = 0.0;
+
+    /*--- Handle averaging and division by volume in one constant. ---*/
+
+    su2double halfOnVol = 0.5 / (nodes->GetVolume(iPoint) + nodes->GetPeriodicVolume(iPoint));
+
+    /*--- Add a contribution due to each neighbor. ---*/
+
+    for (size_t iNeigh = 0; iNeigh < nodes->GetnPoint(iPoint); ++iNeigh) {
+      size_t iEdge = nodes->GetEdge(iPoint,iNeigh);
+      size_t jPoint = nodes->GetPoint(iPoint,iNeigh);
+
+      /*--- Determine if edge points inwards or outwards of iPoint.
+        *    If inwards we need to flip the area vector. ---*/
+
+      su2double dir = (iPoint < jPoint)? 1.0 : -1.0;
+      su2double weight = dir * halfOnVol;
+
+      const auto area = geometry.edges->GetNormal(iEdge);
+
+      for (size_t jDim = 0; jDim < nDim; ++jDim) {
+        for (size_t iVar = varBegin; iVar < varEnd; ++iVar) {
+          su2double flux = weight * (gradient(iPoint, iVar, jDim) + gradient(jPoint, iVar, jDim));
+          for (size_t iDim = 0; iDim < nDim; ++iDim) {
+            size_t ind = (iDim <= jDim) ? iDim*nDim - ((iDim - 1)*iDim)/2 + jDim - iDim
+                                        : jDim*nDim - ((jDim - 1)*jDim)/2 + iDim - jDim;
+            diagScale = (iDim == jDim)? 1.0 : 0.5;
+            hessian(iPoint, iVar, ind) += diagScale * flux * area[iDim];
+          } // iDims
+        } // variables
+      } // jDims
+    } // neighbors
+  } // points
+  END_SU2_OMP_FOR
+
+  /*--- Add edges of markers that contribute to the Hessians ---*/
+  for (size_t iMarker = 0; iMarker < geometry.GetnMarker(); ++iMarker) {
+    if ((config.GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) &&
+      (config.GetMarker_All_KindBC(iMarker) != NEARFIELD_BOUNDARY) &&
+      (config.GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) {
+    /*--- Work is shared in inner loop as two markers
+      *    may try to update the same point. ---*/
+
+      SU2_OMP_FOR_STAT(32)
+      for (size_t iVertex = 0; iVertex < geometry.GetnVertex(iMarker); ++iVertex) {
+        size_t iPoint = geometry.vertex[iMarker][iVertex]->GetNode();
+        auto nodes = geometry.nodes;
+
+        /*--- Halo points do not need to be considered. ---*/
+
+        if (!nodes->GetDomain(iPoint)) continue;
+
+        su2double volume = nodes->GetVolume(iPoint) + nodes->GetPeriodicVolume(iPoint);
+        const auto area = geometry.vertex[iMarker][iVertex]->GetNormal();
+
+        for (size_t jDim = 0; jDim < nDim; ++jDim) {
+          for (size_t iVar = varBegin; iVar < varEnd; iVar++) {
+            su2double flux = gradient(iPoint, iVar, jDim) / volume;
+            for (size_t iDim = 0; iDim < nDim; ++iDim) {
+              size_t ind = (iDim <= jDim) ? iDim * nDim - ((iDim - 1) * iDim)/2 + jDim - iDim
+                                          : jDim * nDim - ((jDim - 1) * jDim)/2 + iDim - jDim;
+              diagScale = (iDim == jDim)? 1.0 : 0.5;
+              hessian(iPoint, iVar, ind) -= diagScale * flux * area[iDim];
+            } // iDims
+          } // variables
+        } // jDims
+      } // vertices
+      END_SU2_OMP_FOR
+    } //found right marker
+  } // iMarkers
+
+  /*--- Compute the corrections for symmetry planes and Euler walls. ---*/
+  /*--- TODO? correctHessiansSymmetry<nDim>() ---*/
+
+  /*--- If no solver was provided we do not communicate ---*/
+
+  if (solver == nullptr) return;
+
+  /*--- Account for periodic contributions. ---*/
+
+  for (size_t iPeriodic = 1; iPeriodic <= config.GetnMarker_Periodic()/2; ++iPeriodic)
+  {
+    solver->InitiatePeriodicComms(&geometry, &config, iPeriodic, kindPeriodicComm);
+    solver->CompletePeriodicComms(&geometry, &config, iPeriodic, kindPeriodicComm);
+  }
+
+  /*--- Obtain the gradients at halo points from the MPI ranks that own them. ---*/
+
+  solver->InitiateComms(&geometry, &config, kindMpiComm);
+  solver->CompleteComms(&geometry, &config, kindMpiComm);
+
+}
 }  // namespace detail
 
 
@@ -205,3 +322,23 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER
       break;
   }
 }
+
+template<class GradientType>
+void computeHessiansGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PERIODIC_QUANTITIES kindPeriodicComm,
+                               CGeometry& geometry, const CConfig& config, const GradientType& gradient,
+                               const size_t varBegin, const size_t varEnd, const int idxVel, GradientType& hessian) {
+  switch (geometry.GetnDim()) {
+  case 2:
+    detail::computeHessiansGreenGauss<2>(solver, kindMpiComm, kindPeriodicComm, geometry, config, gradient,
+                                         varBegin, varEnd, idxVel, hessian);
+    break;
+  case 3:
+    detail::computeHessiansGreenGauss<3>(solver, kindMpiComm, kindPeriodicComm, geometry, config, gradient,
+                                         varBegin, varEnd, idxVel, hessian);
+    break;
+  default:
+    SU2_MPI::Error("Too many dimensions to compute Hessians.", CURRENT_FUNCTION);
+    break;
+  }
+}
+
diff --git a/SU2_CFD/include/gradients/computeMetrics.hpp b/SU2_CFD/include/gradients/computeMetrics.hpp
new file mode 100644
index 000000000000..8917d580973b
--- /dev/null
+++ b/SU2_CFD/include/gradients/computeMetrics.hpp
@@ -0,0 +1,374 @@
+/*!
+ * \file computeMetrics.hpp
+ * \brief Generic implementation of the metric tensor computation.
+ * \note This allows the same implementation to be used for goal-oriented
+ *       or feature-based mesh adaptation.
+ * \author B. Munguía
+ * \version 8.3.0 "Harrier"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md)
+ *
+ * SU2 is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * SU2 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <vector>
+#include <algorithm>
+#include <limits>
+#include <cmath>
+#include "../../../Common/include/parallelization/omp_structure.hpp"
+#include "../../../Common/include/linear_algebra/blas_structure.hpp"
+#include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
+
+namespace tensor {
+  struct metric {
+    template <class MatrixType, class MetricType>
+    static void get(MetricType& metric_field, const unsigned long iPoint,
+                    const unsigned short iSensor, MatrixType& mat, unsigned short nDim) {
+      switch(nDim) {
+        case 2: {
+          mat[0][0] = metric_field(iPoint, 0); mat[0][1] = metric_field(iPoint, 1);
+          mat[1][0] = metric_field(iPoint, 1); mat[1][1] = metric_field(iPoint, 2);
+          break;
+        }
+        case 3: {
+          mat[0][0] = metric_field(iPoint, 0); mat[0][1] = metric_field(iPoint, 1); mat[0][2] = metric_field(iPoint, 2);
+          mat[1][0] = metric_field(iPoint, 1); mat[1][1] = metric_field(iPoint, 3); mat[1][2] = metric_field(iPoint, 4);
+          mat[2][0] = metric_field(iPoint, 2); mat[2][1] = metric_field(iPoint, 4); mat[2][2] = metric_field(iPoint, 5);
+          break;
+        }
+      }
+    }
+
+    template <class ScalarType, class MatrixType, class MetricType>
+    static void set(MetricType& metric_field, const unsigned long iPoint,
+                    const unsigned short iSensor, MatrixType& mat, ScalarType scale,
+                    unsigned short nDim) {
+      switch(nDim) {
+        case 2: {
+          metric_field(iPoint, 0) = mat[0][0] * scale;
+          metric_field(iPoint, 1) = mat[0][1] * scale;
+          metric_field(iPoint, 2) = mat[1][1] * scale;
+          break;
+        }
+        case 3: {
+          metric_field(iPoint, 0) = mat[0][0] * scale;
+          metric_field(iPoint, 1) = mat[0][1] * scale;
+          metric_field(iPoint, 2) = mat[0][2] * scale;
+          metric_field(iPoint, 3) = mat[1][1] * scale;
+          metric_field(iPoint, 4) = mat[1][2] * scale;
+          metric_field(iPoint, 5) = mat[2][2] * scale;
+          break;
+        }
+      }
+    }
+  };
+
+  struct hessian {
+    template <class MatrixType, class MetricType>
+    static void get(MetricType& metric_field, const unsigned long iPoint,
+                    const unsigned short iSensor, MatrixType& mat, unsigned short nDim) {
+      switch(nDim) {
+        case 2: {
+          mat[0][0] = metric_field(iPoint, iSensor, 0); mat[0][1] = metric_field(iPoint, iSensor, 1);
+          mat[1][0] = metric_field(iPoint, iSensor, 1); mat[1][1] = metric_field(iPoint, iSensor, 2);
+          break;
+        }
+        case 3: {
+          mat[0][0] = metric_field(iPoint, iSensor, 0); mat[0][1] = metric_field(iPoint, iSensor, 1); mat[0][2] = metric_field(iPoint, iSensor, 2);
+          mat[1][0] = metric_field(iPoint, iSensor, 1); mat[1][1] = metric_field(iPoint, iSensor, 3); mat[1][2] = metric_field(iPoint, iSensor, 4);
+          mat[2][0] = metric_field(iPoint, iSensor, 2); mat[2][1] = metric_field(iPoint, iSensor, 4); mat[2][2] = metric_field(iPoint, iSensor, 5);
+          break;
+        }
+      }
+    }
+
+    template <class ScalarType, class MatrixType, class MetricType>
+    static void set(MetricType& metric_field, const unsigned long iPoint,
+                    const unsigned short iSensor,  MatrixType& mat, ScalarType scale,
+                    unsigned short nDim) {
+      switch(nDim) {
+        case 2: {
+          metric_field(iPoint, iSensor, 0) = mat[0][0] * scale;
+          metric_field(iPoint, iSensor, 1) = mat[0][1] * scale;
+          metric_field(iPoint, iSensor, 2) = mat[1][1] * scale;
+          break;
+        }
+        case 3: {
+          metric_field(iPoint, iSensor, 0) = mat[0][0] * scale;
+          metric_field(iPoint, iSensor, 1) = mat[0][1] * scale;
+          metric_field(iPoint, iSensor, 2) = mat[0][2] * scale;
+          metric_field(iPoint, iSensor, 3) = mat[1][1] * scale;
+          metric_field(iPoint, iSensor, 4) = mat[1][2] * scale;
+          metric_field(iPoint, iSensor, 5) = mat[2][2] * scale;
+          break;
+        }
+      }
+    }
+  };
+}
+
+namespace detail {
+
+/*!
+ * \brief Compute determinant of eigenvalues for different dimensions.
+ * \param[in] EigVal - Array of eigenvalues.
+ * \return Determinant value.
+ */
+template<size_t nDim, class ScalarType>
+ScalarType computeDeterminant(const ScalarType* EigVal) {
+  if constexpr (nDim == 2) {
+    return EigVal[0] * EigVal[1];
+  } else if constexpr (nDim == 3) {
+    return EigVal[0] * EigVal[1] * EigVal[2];
+  } else {
+    static_assert(nDim == 2 || nDim == 3, "Only 2D and 3D supported");
+    return ScalarType(0.0);
+  }
+}
+
+/*!
+ * \brief Make the eigenvalues of the metrics positive.
+ * \param[in] geometry - Geometrical definition of the problem.
+ * \param[in] config - Definition of the particular problem.
+ * \param[in] iSensor - Index of the sensor to work on.
+ * \param[in,out] metric - Metric container.
+ */
+template<size_t nDim, class ScalarType, class Tensor, class MetricType>
+void setPositiveDefiniteMetrics(CGeometry& geometry, const CConfig& config,
+                               unsigned short iSensor, MetricType& metric) {
+
+  const unsigned long nPointDomain = geometry.GetnPointDomain();
+
+  ScalarType A[nDim][nDim], EigVec[nDim][nDim], EigVal[nDim], work[nDim];
+
+  /*--- Minimum eigenvalue threshold ---*/
+  const ScalarType eps = 1e-12;
+
+  for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) {
+    /*--- Get full metric tensor ---*/
+    Tensor::get(metric, iPoint, iSensor, A, nDim);
+
+    /*--- Compute eigenvalues and eigenvectors ---*/
+    CBlasStructure::EigenDecomposition(A, EigVec, EigVal, nDim, work);
+
+    /*--- Make positive definite by taking absolute value of eigenvalues ---*/
+    /*--- Handle NaN and very small values that could cause numerical issues ---*/
+    for (auto iDim = 0; iDim < nDim; iDim++) {
+      if (std::isnan(EigVal[iDim])) {
+        /*--- NaN detected, set to small positive value ---*/
+        EigVal[iDim] = eps;
+      } else {
+        /*--- Take absolute value and ensure minimum threshold ---*/
+        EigVal[iDim] = max(fabs(EigVal[iDim]), eps);
+      }
+    }
+
+    CBlasStructure::EigenRecomposition(A, EigVec, EigVal, nDim);
+
+    /*--- Store upper half of metric tensor ---*/
+    Tensor::set(metric, iPoint, iSensor, A, 1.0, nDim);
+  }
+}
+
+/*!
+ * \brief Integrate the Hessian field for the Lp-norm normalization of the metric.
+ * \param[in] geometry - Geometrical definition of the problem.
+ * \param[in] config - Definition of the particular problem.
+ * \param[in] iSensor - Index of the sensor to work on.
+ * \param[in] metric - Metric container.
+ * \return Integral of the metric tensor determinant.
+*/
+template<size_t nDim, class ScalarType, class MetricType>
+ScalarType integrateMetrics(CGeometry& geometry, const CConfig& config,
+                            unsigned short iSensor, MetricType& metric) {
+
+  const unsigned long nPointDomain = geometry.GetnPointDomain();
+
+  /*--- Constants defining normalization ---*/
+  const ScalarType p = config.GetMetric_Norm();
+  const ScalarType normExp = p / (2.0 * p + nDim);
+
+  ScalarType localIntegral = 0.0;
+  ScalarType globalIntegral = 0.0;
+  for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) {
+    auto nodes = geometry.nodes;
+
+    /*--- Calculate determinant ---*/
+    ScalarType det;
+    if constexpr (nDim == 2) {
+      const ScalarType m00 = metric(iPoint, 0);
+      const ScalarType m01 = metric(iPoint, 1);
+      const ScalarType m11 = metric(iPoint, 2);
+      det = m00 * m11 - m01 * m01;
+    } else if constexpr (nDim == 3) {
+      const ScalarType m00 = metric(iPoint, 0);
+      const ScalarType m01 = metric(iPoint, 1);
+      const ScalarType m02 = metric(iPoint, 2);
+      const ScalarType m11 = metric(iPoint, 3);
+      const ScalarType m12 = metric(iPoint, 4);
+      const ScalarType m22 = metric(iPoint, 5);
+      det = m00 * (m11 * m22 - m12 * m12) - m01 * (m01 * m22 - m02 * m12) + m02 * (m01 * m12 - m02 * m11);
+    }
+
+    /*--- Integrate determinant ---*/
+    const ScalarType Vol = SU2_TYPE::GetValue(nodes->GetVolume(iPoint));
+    localIntegral += pow(abs(det), normExp) * Vol;
+  }
+
+  CBaseMPIWrapper::Allreduce(&localIntegral, &globalIntegral, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+
+  return globalIntegral;
+}
+
+/*!
+ * \brief Perform an Lp-norm normalization of the metric.
+ * \param[in] geometry - Geometrical definition of the problem.
+ * \param[in] config - Definition of the particular problem.
+ * \param[in] iSensor - Index of the sensor to work on.
+ * \param[in] integral - Integral of the metric tensor determinant.
+ * \param[in,out] metric - Metric container.
+ */
+template<size_t nDim, class ScalarType, class Tensor, class MetricType>
+void normalizeMetrics(CGeometry& geometry, const CConfig& config,
+                      unsigned short iSensor, ScalarType integral,
+                      MetricType& metric) {
+
+  const unsigned long nPointDomain = geometry.GetnPointDomain();
+
+  /*--- Constants defining normalization ---*/
+  const ScalarType p = config.GetMetric_Norm();
+  const ScalarType N = SU2_TYPE::GetValue(config.GetMetric_Complexity());
+  const ScalarType globalFactor = pow(N / integral, 2.0 / nDim);
+  const ScalarType normExp = -1.0 / (2.0 * p + nDim);
+
+  /*--- Size constraints ---*/
+  const ScalarType hmin = SU2_TYPE::GetValue(config.GetMetric_Hmin());
+  const ScalarType hmax = SU2_TYPE::GetValue(config.GetMetric_Hmax());
+  const ScalarType eigmax = 1.0 / pow(hmin, 2.0);
+  const ScalarType eigmin = 1.0 / pow(hmax, 2.0);
+  const ScalarType armax2 = pow(SU2_TYPE::GetValue(config.GetMetric_ARmax()), 2.0);
+
+  ScalarType A[nDim][nDim], EigVec[nDim][nDim], EigVal[nDim], work[nDim];
+
+  for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) {
+    /*--- Decompose metric ---*/
+    Tensor::get(metric, iPoint, iSensor, A, nDim);
+    CBlasStructure::EigenDecomposition(A, EigVec, EigVal, nDim, work);
+
+    /*--- Normalize eigenvalues ---*/
+    const ScalarType det = computeDeterminant<nDim>(EigVal);
+    const ScalarType factor = globalFactor * pow(abs(det), normExp);
+    for (auto iDim = 0u; iDim < nDim; ++iDim)
+      EigVal[iDim] = factor * EigVal[iDim];
+
+    /*--- Clip by user-specified size constraints ---*/
+    for (auto iDim = 0u; iDim < nDim; ++iDim)
+      EigVal[iDim] = min(max(abs(EigVal[iDim]), eigmin), eigmax);
+
+    /*--- Clip by user-specified aspect ratio ---*/
+    unsigned short iMax = 0;
+    for (auto iDim = 1; iDim < nDim; ++iDim)
+      iMax = (EigVal[iDim] > EigVal[iMax])? iDim : iMax;
+
+    for (auto iDim = 0u; iDim < nDim; ++iDim)
+      EigVal[iDim] = max(EigVal[iDim], EigVal[iMax]/armax2);
+
+    /*--- Recompose and store metric ---*/
+    CBlasStructure::EigenRecomposition(A, EigVec, EigVal, nDim);
+    Tensor::set(metric, iPoint, iSensor, A, 1.0, nDim);
+  }
+}
+} // end namespace detail
+
+/*!
+ * \brief Make the eigenvalues of the metrics positive.
+ * \param[in] geometry - Geometrical definition of the problem.
+ * \param[in] config - Definition of the particular problem.
+ * \param[in] iSensor - Index of the sensor to work on.
+ * \param[in,out] metric - Metric container.
+ */
+template<class ScalarType, class Tensor, class MetricType>
+void setPositiveDefiniteMetrics(CGeometry& geometry, const CConfig& config,
+                               unsigned short iSensor, MetricType& metric) {
+  switch (geometry.GetnDim()) {
+    case 2:
+      detail::setPositiveDefiniteMetrics<2, ScalarType, Tensor>(geometry, config, iSensor, metric);
+      break;
+    case 3:
+      detail::setPositiveDefiniteMetrics<3, ScalarType, Tensor>(geometry, config, iSensor, metric);
+      break;
+    default:
+      SU2_MPI::Error("Too many dimensions for metric computation.", CURRENT_FUNCTION);
+      break;
+  }
+}
+
+/*!
+ * \brief Integrate the Hessian field for the Lp-norm normalization of the metric.
+ * \param[in] geometry - Geometrical definition of the problem.
+ * \param[in] config - Definition of the particular problem.
+ * \param[in] iSensor - Index of the sensor to work on.
+ * \param[in] metric - Metric container.
+ * \return Integral of the metric tensor determinant.
+*/
+template<class ScalarType, class MetricType>
+ScalarType integrateMetrics(CGeometry& geometry, const CConfig& config,
+                            unsigned short iSensor, MetricType& metric) {
+  su2double integral;
+  switch (geometry.GetnDim()) {
+    case 2:
+      integral = detail::integrateMetrics<2, ScalarType>(geometry, config, iSensor, metric);
+      break;
+    case 3:
+      integral = detail::integrateMetrics<3, ScalarType>(geometry, config, iSensor, metric);
+      break;
+    default:
+      SU2_MPI::Error("Too many dimensions for metric integration.", CURRENT_FUNCTION);
+      break;
+  }
+
+  return integral;
+}
+
+/*!
+ * \brief Perform an Lp-norm normalization of the metric.
+ * \param[in] geometry - Geometrical definition of the problem.
+ * \param[in] config - Definition of the particular problem.
+ * \param[in] iSensor - Index of the sensor to work on.
+ * \param[in] integral - Integral of the metric tensor determinant.
+ * \param[in,out] metric - Metric container.
+ */
+template<class ScalarType, class Tensor, class MetricType>
+void normalizeMetrics(CGeometry& geometry, const CConfig& config,
+                     unsigned short iSensor, ScalarType integral,
+                     MetricType& metric) {
+  switch (geometry.GetnDim()) {
+    case 2:
+      detail::normalizeMetrics<2, ScalarType, Tensor>(geometry, config, iSensor, integral, metric);
+      break;
+    case 3:
+      detail::normalizeMetrics<3, ScalarType, Tensor>(geometry, config, iSensor, integral, metric);
+      break;
+    default:
+      SU2_MPI::Error("Too many dimensions for metric normalization.", CURRENT_FUNCTION);
+      break;
+  }
+}
diff --git a/SU2_CFD/include/output/CFlowOutput.hpp b/SU2_CFD/include/output/CFlowOutput.hpp
index 63867b091ae4..7c8c23af181f 100644
--- a/SU2_CFD/include/output/CFlowOutput.hpp
+++ b/SU2_CFD/include/output/CFlowOutput.hpp
@@ -352,4 +352,20 @@ class CFlowOutput : public CFVMOutput{
    * \param[in] config - Definition of the particular problem per zone.
    */
   void SetFixedCLScreenOutput(const CConfig *config);
+
+  /*!
+   * \brief Add mesh adaptation outputs.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void AddMeshAdaptationOutputs(const CConfig* config);
+
+  /*!
+   * \brief Load mesh adaptation outputs.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] solver - The container holding all solution data.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] iPoint - Index of the point.
+   */
+  void LoadMeshAdaptationOutputs(const CConfig* config, const CSolver* const* solver, const CGeometry* geometry,
+                                 unsigned long iPoint);
 };
diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp
index e3f029a37538..d8a4286d20cc 100644
--- a/SU2_CFD/include/solvers/CEulerSolver.hpp
+++ b/SU2_CFD/include/solvers/CEulerSolver.hpp
@@ -273,6 +273,57 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMP
   virtual unsigned long SetPrimitive_Variables(CSolver **solver_container,
                                                const CConfig *config);
 
+  /*!
+   * \brief Store the primitive variables needed for adaptation.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetPrimitive_Adapt(CGeometry *geometry, const CConfig *config) final {
+    const auto nSensor = config->GetnMetric_Sensor();
+    for (auto iSensor = 0; iSensor < nSensor; iSensor++) {
+      switch (config->GetMetric_Sensor(iSensor)) {
+        case METRIC_SENSOR::DENSITY:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = nodes->GetDensity(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+        case METRIC_SENSOR::PRESSURE:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = nodes->GetPressure(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+        case METRIC_SENSOR::TOTAL_PRESSURE:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double Mach = sqrt(nodes->GetVelocity2(iPoint))/nodes->GetSoundSpeed(iPoint);
+            const su2double prim_var = nodes->GetPressure(iPoint)*pow((1+0.2*Mach*Mach), 3.5);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+        case METRIC_SENSOR::TEMPERATURE:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = nodes->GetTemperature(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+        case METRIC_SENSOR::ENERGY:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = nodes->GetEnergy(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+        case METRIC_SENSOR::MACH:
+        default:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = sqrt(nodes->GetVelocity2(iPoint)) / nodes->GetSoundSpeed(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+      }
+    }
+  }
+
   /*!
    * \brief Set gradients of coefficients for fixed CL mode
    * \param[in] config - Definition of the particular problem.
diff --git a/SU2_CFD/include/solvers/CIncEulerSolver.hpp b/SU2_CFD/include/solvers/CIncEulerSolver.hpp
index 54fd616d0de8..941f46c2eb52 100644
--- a/SU2_CFD/include/solvers/CIncEulerSolver.hpp
+++ b/SU2_CFD/include/solvers/CIncEulerSolver.hpp
@@ -86,6 +86,50 @@ class CIncEulerSolver : public CFVMFlowSolverBase<CIncEulerVariable, ENUM_REGIME
    */
   virtual unsigned long SetPrimitive_Variables(CSolver **solver_container, const CConfig *config);
 
+  /*!
+   * \brief Store the primitive variables needed for adaptation.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetPrimitive_Adapt(CGeometry *geometry, const CConfig *config) final {
+    const auto nSensor = config->GetnMetric_Sensor();
+    for (auto iSensor = 0; iSensor < nSensor; iSensor++) {
+      switch (config->GetMetric_Sensor(iSensor)) {
+        case METRIC_SENSOR::DENSITY:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = nodes->GetDensity(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+        case METRIC_SENSOR::TOTAL_PRESSURE:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = nodes->GetPressure(iPoint) + 0.5*nodes->GetDensity(iPoint)*nodes->GetVelocity2(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+        case METRIC_SENSOR::TEMPERATURE:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = nodes->GetTemperature(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+        case METRIC_SENSOR::ENERGY:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = nodes->GetEnergy(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+        case METRIC_SENSOR::PRESSURE:
+        default:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = nodes->GetPressure(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+      }
+    }
+  }
+
   /*!
    * \brief Update the Beta parameter for the incompressible preconditioner.
    * \param[in] geometry - Geometrical definition of the problem.
diff --git a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp
index 9077ea195d89..3b0dc0d8389a 100644
--- a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp
+++ b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp
@@ -195,6 +195,56 @@ class CNEMOEulerSolver : public CFVMFlowSolverBase<CNEMOEulerVariable, ENUM_REGI
   virtual unsigned long SetPrimitive_Variables(CSolver **solver_container,
                                                CConfig *config, bool Output);
 
+  /*!
+   * \brief Store the primitive variables needed for adaptation.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetPrimitive_Adapt(CGeometry *geometry, const CConfig *config) final {
+    const auto nSensor = config->GetnMetric_Sensor();
+    for (auto iSensor = 0; iSensor < nSensor; iSensor++) {
+      switch (config->GetMetric_Sensor(iSensor)) {
+        case METRIC_SENSOR::PRESSURE:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = nodes->GetPressure(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+        case METRIC_SENSOR::TEMPERATURE:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = nodes->GetTemperature(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+        case METRIC_SENSOR::TEMPERATURE_VE:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = nodes->GetTemperature_ve(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+        case METRIC_SENSOR::ENERGY:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = nodes->GetEnergy(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+        case METRIC_SENSOR::ENERGY_VE:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = *nodes->GetEve(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+        case METRIC_SENSOR::MACH:
+        default:
+          for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+            const su2double prim_var = sqrt(nodes->GetVelocity2(iPoint)) / nodes->GetSoundSpeed(iPoint);
+            nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var);
+          }
+          break;
+      }
+    }
+  }
+
   /*!
    * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over
    * a nonlinear iteration for stability.
diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp
index f0b3a4c493d1..9f801d22925a 100644
--- a/SU2_CFD/include/solvers/CSolver.hpp
+++ b/SU2_CFD/include/solvers/CSolver.hpp
@@ -85,7 +85,8 @@ class CSolver {
   nSecondaryVar,                 /*!< \brief Number of primitive variables of the problem. */
   nSecondaryVarGrad,             /*!< \brief Number of primitive variables of the problem in the gradient computation. */
   nVarGrad,                      /*!< \brief Number of variables for deallocating the LS Cvector. */
-  nDim;                          /*!< \brief Number of dimensions of the problem. */
+  nDim,                          /*!< \brief Number of dimensions of the problem. */
+  nSymMat;                       /*!< \brief Number of symmetric matrix componenents for Hessian and metric tensor. */
   unsigned long nPoint;          /*!< \brief Number of points of the computational grid. */
   unsigned long nPointDomain;    /*!< \brief Number of points of the computational grid. */
   su2double Max_Delta_Time, /*!< \brief Maximum value of the delta time for all the control volumes. */
@@ -580,6 +581,22 @@ class CSolver {
    */
   inline virtual void SetPrimitive_Limiter(CGeometry *geometry, const CConfig *config) { }
 
+   /*!
+   * \brief A virtual member.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  virtual void SetPrimitive_Adapt(CGeometry *geometry, const CConfig *config) { }
+
+  /*!
+   * \brief Compute the Green-Gauss Hessian of the solution.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] idxVel - Index to velocity, -1 if no velocity is present in the solver.
+   * \param[in] reconstruction - indicator that the gradient being computed is for upwind reconstruction.
+   */
+  void SetHessian_GG(CGeometry *geometry, const CConfig *config, short idxVel, const unsigned short Kind_Solver);
+
   /*!
    * \brief Compute the projection of a variable for MUSCL reconstruction.
    * \note The result should be halved when added to i (or subtracted from j).
@@ -4349,24 +4366,42 @@ class CSolver {
     END_SU2_OMP_FOR
   }
 
-inline void CustomSourceResidual(CGeometry *geometry, CSolver **solver_container,
-                                 CNumerics **numerics_container, CConfig *config, unsigned short iMesh) {
+  inline void CustomSourceResidual(CGeometry *geometry, CSolver **solver_container,
+                                  CNumerics **numerics_container, CConfig *config, unsigned short iMesh) {
 
-  AD::StartNoSharedReading();
+    AD::StartNoSharedReading();
 
-  SU2_OMP_FOR_STAT(roundUpDiv(nPointDomain,2*omp_get_max_threads()))
-  for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) {
-    /*--- Get control volume size. ---*/
-    su2double Volume = geometry->nodes->GetVolume(iPoint);
-    /*--- Compute the residual for this control volume and subtract. ---*/
-    for (auto iVar = 0ul; iVar < nVar; iVar++) {
-      LinSysRes(iPoint,iVar) -= base_nodes->GetUserDefinedSource()(iPoint, iVar) * Volume;
+    SU2_OMP_FOR_STAT(roundUpDiv(nPointDomain,2*omp_get_max_threads()))
+    for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) {
+      /*--- Get control volume size. ---*/
+      su2double Volume = geometry->nodes->GetVolume(iPoint);
+      /*--- Compute the residual for this control volume and subtract. ---*/
+      for (auto iVar = 0ul; iVar < nVar; iVar++) {
+        LinSysRes(iPoint,iVar) -= base_nodes->GetUserDefinedSource()(iPoint, iVar) * Volume;
+      }
     }
+    END_SU2_OMP_FOR
+
+    AD::EndNoSharedReading();
   }
-  END_SU2_OMP_FOR
 
-  AD::EndNoSharedReading();
-}
+  /*!
+   * \brief Compute the metric tensor field.
+   * \param[in] solver - Physical definition of the problem.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void ComputeMetric(CSolver **solver, CGeometry *geometry, const CConfig *config);
+
+  /*!
+   * \brief Add contribution to the metric field.
+   * \param[in] solver - Physical definition of the problem.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iSensor - Index of the sensor to work on.
+   */
+  void AddMetrics(CSolver **solver, const CGeometry *geometry, const CConfig *config,
+                  const unsigned short iSensor);
 
 protected:
   /*!
diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp
index dc905d79fcb4..5cd23008e271 100644
--- a/SU2_CFD/include/variables/CVariable.hpp
+++ b/SU2_CFD/include/variables/CVariable.hpp
@@ -103,6 +103,11 @@ class CVariable {
 
   VectorType SolutionExtra_BGS_k; /*!< \brief Intermediate storage, enables cross term extraction as that is also pushed to Solution. */
 
+  MatrixType Primitive_Adapt;     /*!< \brief Variables for which we need gradients for anisotropy in mesh adaptation. */
+  CVectorOfMatrix Gradient_Adapt; /*!< \brief Gradient of sensor used for anisotropy in mesh adaptation. */
+  CVectorOfMatrix Hessian;        /*!< \brief Hessian of sensor used for anisotropy in mesh adaptation. */
+  su2matrix<double> Metric;       /*!< \brief Metric tensor used for anisotropy in mesh adaptation. */
+
  protected:
   unsigned long nPoint = 0;  /*!< \brief Number of points in the domain. */
   unsigned long nDim = 0;      /*!< \brief Number of dimension of the problem. */
@@ -112,6 +117,8 @@ class CVariable {
   unsigned long nSecondaryVar = 0;     /*!< \brief Number of secondary variables. */
   unsigned long nSecondaryVarGrad = 0;   /*!< \brief Number of secondaries for which a gradient is computed. */
   unsigned long nAuxVar = 0; /*!< \brief Number of auxiliary variables. */
+  unsigned long nSymMat = 0;        /*!< \brief Number of symmetric matrix componenents for Hessian and metric tensor. */
+
 
   /*--- Only allow default construction by derived classes. ---*/
   CVariable() = default;
@@ -2371,4 +2378,115 @@ class CVariable {
 
   inline virtual const su2double *GetScalarSources(unsigned long iPoint) const { return nullptr; }
   inline virtual const su2double *GetScalarLookups(unsigned long iPoint) const { return nullptr; }
+
+  /*!
+   * \brief Set the gradient of the solution.
+   * \param[in] iPoint - Point index.
+   * \param[in] gradient - Gradient of the solution.
+   */
+  inline void SetPrimitive_Adapt(unsigned long iPoint, unsigned long iVar, su2double primitive) { Primitive_Adapt(iPoint,iVar) = primitive; }
+
+  /*!
+   * \brief Get the gradient of the entire solution.
+   * \return Reference to gradient.
+   */
+  inline const MatrixType& GetPrimitive_Adapt(void) const { return Primitive_Adapt; }
+
+  /*!
+   * \brief Get the value of the solution gradient.
+   * \param[in] iPoint - Point index.
+   * \return Value of the gradient solution.
+   */
+  inline su2double *GetPrimitive_Adapt(unsigned long iPoint) { return Primitive_Adapt[iPoint]; }
+
+  /*!
+   * \brief Get the value of the solution gradient.
+   * \param[in] iPoint - Point index.
+   * \param[in] iVar - Variable index.
+   * \return Value of the solution gradient.
+   */
+  inline su2double GetPrimitive_Adapt(unsigned long iPoint, unsigned long iVar) const { return Primitive_Adapt(iPoint,iVar); }
+
+  /*!
+   * \brief Set the gradient of the solution.
+   * \param[in] iPoint - Point index.
+   * \param[in] gradient - Gradient of the solution.
+   */
+  inline void SetGradient_Adapt(unsigned long iPoint, su2double** gradient) {
+    for (unsigned long iVar = 0; iVar < nVar; iVar++)
+      for (unsigned long iDim = 0; iDim < nDim; iDim++)
+        Gradient_Adapt(iPoint,iVar,iDim) = gradient[iVar][iDim];
+  }
+
+  /*!
+   * \brief Get the gradient of the entire solution.
+   * \return Reference to gradient.
+   */
+  inline CVectorOfMatrix& GetGradient_Adapt(void) { return Gradient_Adapt; }
+
+  /*!
+   * \brief Get the value of the solution gradient.
+   * \param[in] iPoint - Point index.
+   * \return Value of the gradient solution.
+   */
+  inline CMatrixView<su2double> GetGradient_Adapt(unsigned long iPoint) { return Gradient_Adapt[iPoint]; }
+
+  /*!
+   * \brief Get the value of the solution gradient.
+   * \param[in] iPoint - Point index.
+   * \param[in] iVar - Variable index.
+   * \param[in] iDim - Dimension index.
+   * \return Value of the solution gradient.
+   */
+  inline su2double GetGradient_Adapt(unsigned long iPoint, unsigned long iVar, unsigned long iDim) const { return Gradient_Adapt(iPoint,iVar,iDim); }
+
+  /*!
+   * \brief Set the Hessian of the solution.
+   * \param[in] iPoint - Point index.
+   * \param[in] iVar - Variable index.
+   * \param[in] iMat - Hessian tensor index.
+   * \param[in] hess - Hessian of the solution.
+   */
+  inline void SetHessian(unsigned long iPoint, unsigned long iVar, unsigned long iMat, su2double hess) { Hessian(iPoint,iVar,iMat) = hess; }
+
+  /*!
+   * \brief Get the Hessian tensor field.
+   * \return Reference to the Hessian field.
+   */
+  inline CVectorOfMatrix& GetHessian(void) { return Hessian; }
+
+  /*!
+   * \brief Get the value of the Hessian.
+   * \param[in] iPoint - Point index.
+   * \param[in] iVar - Variable index.
+   * \param[in] iMat - Hessian tensor index.
+   * \return Value of the Hessian.
+   */
+  inline su2double GetHessian(unsigned long iPoint, unsigned long iVar, unsigned long iMat) const { return Hessian(iPoint,iVar,iMat); }
+
+  /*!
+   * \brief Set the value of the metric.
+   * \param[in] iMat - Metric tensor index.
+   * \param[in] metric - Metric value.
+   */
+  inline void SetMetric(unsigned long iPoint, unsigned short iMat, double metric) { Metric(iPoint,iMat) = metric; }
+
+  /*!
+   * \brief Get the metric of the entire solution.
+   * \return Reference to the metric tensor field.
+   */
+  inline su2matrix<double>& GetMetric(void) { return Metric; }
+
+  /*!
+   * \brief Add the value of the metric.
+   * \param[in] iMat - Metric tensor index.
+   * \param[in] metric - Metric value.
+   */
+  inline void AddMetric(unsigned long iPoint, unsigned short iMat, double metric) { Metric(iPoint,iMat) += metric; }
+
+  /*!
+   * \brief Get the value of the metric.
+   * \param[in] iMat  - Metric tensor index.
+   */
+  inline double GetMetric(unsigned long iPoint, unsigned short iMat) const { return Metric(iPoint,iMat); }
 };
diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp
index 1270eef237ab..f4343ba9d150 100644
--- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp
+++ b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp
@@ -175,6 +175,11 @@ void CSinglezoneDriver::Postprocess() {
     iteration_container[ZONE_0][INST_0]->Relaxation(output_container[ZONE_0], integration_container, geometry_container, solver_container,
         numerics_container, config_container, surface_movement, grid_movement, FFDBox, ZONE_0, INST_0);
 
+  /*--- Compute metric for anisotropic mesh adaptation ---*/
+
+  if (config_container[ZONE_0]->GetCompute_Metric())
+    ComputeMetricField();
+
 }
 
 void CSinglezoneDriver::Update() {
@@ -312,3 +317,33 @@ bool CSinglezoneDriver::Monitor(unsigned long TimeIter){
 bool CSinglezoneDriver::GetTimeConvergence() const{
   return output_container[ZONE_0]->GetCauchyCorrectedTimeConvergence(config_container[ZONE_0]);
 }
+
+void CSinglezoneDriver::ComputeMetricField() {
+
+  auto solver = solver_container[ZONE_0][INST_0][MESH_0];
+  auto solver_flow = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL];
+  auto geometry = geometry_container[ZONE_0][INST_0][MESH_0];
+  auto config = config_container[ZONE_0];
+  int idxVel = -1;
+
+  if (rank == MASTER_NODE){
+    cout << endl <<"----------------------------- Compute Metric ----------------------------" << endl;
+    cout << "Storing primitive variables needed for gradients in metric." << endl;
+  }
+  solver_flow->InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION);
+  solver_flow->CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION);
+  solver_flow->Preprocessing(geometry, solver, config, MESH_0, NO_RK_ITER,
+                             RUNTIME_FLOW_SYS, true);
+  solver_flow->SetPrimitive_Adapt(geometry, config);
+
+  if (config->GetKind_Hessian_Method() == GREEN_GAUSS) {
+    if(rank == MASTER_NODE) cout << "Computing Hessians using Green-Gauss." << endl;
+    solver_flow->SetHessian_GG(geometry, config, idxVel, RUNTIME_FLOW_SYS);
+  }
+  else {
+    SU2_MPI::Error("Unsupported Hessian method.", CURRENT_FUNCTION);
+  }
+
+  if(rank == MASTER_NODE) cout << "Computing feature-based metric tensor." << endl;
+  solver_flow->ComputeMetric(solver, geometry, config);
+}
diff --git a/SU2_CFD/src/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp
index 2476498e84d7..6bcae5a287f0 100644
--- a/SU2_CFD/src/output/CFlowCompOutput.cpp
+++ b/SU2_CFD/src/output/CFlowCompOutput.cpp
@@ -303,6 +303,9 @@ void CFlowCompOutput::SetVolumeOutputFields(CConfig *config){
 
   AddCommonFVMOutputs(config);
 
+  // Metric tensor and dual-cell volume
+  AddMeshAdaptationOutputs(config);
+
   if (config->GetTime_Domain()) {
     SetTimeAveragedFields();
   }
diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp
index b355c0698afb..f6c7c38735d4 100644
--- a/SU2_CFD/src/output/CFlowOutput.cpp
+++ b/SU2_CFD/src/output/CFlowOutput.cpp
@@ -4109,3 +4109,41 @@ void CFlowOutput::AddTurboOutput(unsigned short nZone){
   AddHistoryOutput("KineticEnergyLoss_Stage", "KELC_all", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine Kinetic Energy Loss Coefficient", HistoryFieldType::DEFAULT);
   AddHistoryOutput("TotPressureLoss_Stage", "TPLC_all", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine Pressure Loss Coefficient", HistoryFieldType::DEFAULT);
 }
+
+void CFlowOutput::AddMeshAdaptationOutputs(const CConfig* config) {
+
+  // Anisotropic metric tensor, and dual-cell volume
+  if(config->GetCompute_Metric()) {
+    // Common metric components for both 2D and 3D
+    AddVolumeOutput("METRIC_XX", "Metric_xx", "MESH_ADAPT", "x-x-component of the metric");
+    AddVolumeOutput("METRIC_XY", "Metric_xy", "MESH_ADAPT", "x-y-component of the metric");
+    AddVolumeOutput("METRIC_YY", "Metric_yy", "MESH_ADAPT", "y-y-component of the metric");
+
+    // Additional components for 3D
+    if (nDim == 3) {
+      AddVolumeOutput("METRIC_XZ", "Metric_xz", "MESH_ADAPT", "x-z-component of the metric");
+      AddVolumeOutput("METRIC_YZ", "Metric_yz", "MESH_ADAPT", "y-z-component of the metric");
+      AddVolumeOutput("METRIC_ZZ", "Metric_zz", "MESH_ADAPT", "z-z-component of the metric");
+    }
+  }
+}
+
+void CFlowOutput::LoadMeshAdaptationOutputs(const CConfig* config, const CSolver* const* solver, const CGeometry* geometry,
+                                            unsigned long iPoint) {
+
+  const auto* Node_Flow = solver[FLOW_SOL]->GetNodes();
+  if(config->GetCompute_Metric()) {
+    // Common metric components for both 2D and 3D
+    SetVolumeOutputValue("METRIC_XX", iPoint, Node_Flow->GetMetric(iPoint, 0));
+    SetVolumeOutputValue("METRIC_XY", iPoint, Node_Flow->GetMetric(iPoint, 1));
+    SetVolumeOutputValue("METRIC_YY", iPoint, Node_Flow->GetMetric(iPoint, 2));
+
+    // Additional components for 3D
+    if (nDim == 3) {
+      SetVolumeOutputValue("METRIC_XZ", iPoint, Node_Flow->GetMetric(iPoint, 3));
+      SetVolumeOutputValue("METRIC_YZ", iPoint, Node_Flow->GetMetric(iPoint, 4));
+      SetVolumeOutputValue("METRIC_ZZ", iPoint, Node_Flow->GetMetric(iPoint, 5));
+    }
+  }
+}
+
diff --git a/SU2_CFD/src/output/filewriter/CParaviewXMLFileWriter.cpp b/SU2_CFD/src/output/filewriter/CParaviewXMLFileWriter.cpp
index 5a848b2e11a9..da99b870fe4b 100644
--- a/SU2_CFD/src/output/filewriter/CParaviewXMLFileWriter.cpp
+++ b/SU2_CFD/src/output/filewriter/CParaviewXMLFileWriter.cpp
@@ -53,8 +53,10 @@ void CParaviewXMLFileWriter::WriteData(string val_filename){
   }
 
   /*--- We always have 3 coords, independent of the actual value of nDim ---*/
+  /*--- We also always have 6 tensor components ---*/
 
   const int NCOORDS = 3;
+  const int NTENSOR = 6;
   const unsigned short nDim = dataSorter->GetnDim();
   unsigned short iDim = 0;
 
@@ -139,30 +141,85 @@ void CParaviewXMLFileWriter::WriteData(string val_filename){
     fieldname.erase(remove(fieldname.begin(), fieldname.end(), '"'),
                     fieldname.end());
 
-    /*--- Check whether this field is a vector or scalar. ---*/
+    /*--- Check whether this field is a vector, tensor, or scalar. ---*/
+    /*--- Check tensor components first (longer patterns) to avoid false matches ---*/
 
-    bool output_variable = true, isVector = false;
-    size_t found = fieldNames[iField].find("_x");
+    bool output_variable = true, isVector = false, isTensor = false;
+
+    // Check for tensor components first (_xx, _xy, _yy, _xz, _yz, _zz)
+    size_t found = fieldNames[iField].find("_xx");
     if (found!=string::npos) {
       output_variable = true;
-      isVector        = true;
+      isTensor        = true;
+    }
+    found = fieldNames[iField].find("_xy");
+    if (found!=string::npos) {
+      /*--- We have found a tensor, so skip the XY component. ---*/
+      output_variable = false;
+      isTensor        = true;
+      VarCounter++;
+    }
+    found = fieldNames[iField].find("_yy");
+    if (found!=string::npos) {
+      /*--- We have found a tensor, so skip the YY component. ---*/
+      output_variable = false;
+      isTensor        = true;
+      VarCounter++;
+    }
+    found = fieldNames[iField].find("_xz");
+    if (found!=string::npos) {
+      /*--- We have found a tensor, so skip the XZ component. ---*/
+      output_variable = false;
+      isTensor        = true;
+      VarCounter++;
     }
-    found = fieldNames[iField].find("_y");
+    found = fieldNames[iField].find("_yz");
     if (found!=string::npos) {
-      /*--- We have found a vector, so skip the Y component. ---*/
+      /*--- We have found a tensor, so skip the YZ component. ---*/
       output_variable = false;
+      isTensor        = true;
       VarCounter++;
     }
-    found = fieldNames[iField].find("_z");
+    found = fieldNames[iField].find("_zz");
     if (found!=string::npos) {
-      /*--- We have found a vector, so skip the Z component. ---*/
+      /*--- We have found a tensor, so skip the ZZ component. ---*/
       output_variable = false;
+      isTensor        = true;
       VarCounter++;
     }
 
-    /*--- Write the point data as an <X,Y,Z> vector or a scalar. ---*/
+    // Check for vector components only if not a tensor (_x, _y, _z)
+    if (!isTensor) {
+      found = fieldNames[iField].find("_x");
+      if (found!=string::npos) {
+        output_variable = true;
+        isVector        = true;
+      }
+      found = fieldNames[iField].find("_y");
+      if (found!=string::npos) {
+        /*--- We have found a vector, so skip the Y component. ---*/
+        output_variable = false;
+        VarCounter++;
+      }
+      found = fieldNames[iField].find("_z");
+      if (found!=string::npos) {
+        /*--- We have found a vector, so skip the Z component. ---*/
+        output_variable = false;
+        VarCounter++;
+      }
+    }
+
+    /*--- Write the point data as an <X,Y,Z> vector, tensor, or a scalar. ---*/
+
+    if (output_variable && isTensor) {
 
-    if (output_variable && isVector) {
+      /*--- Adjust the string name to remove the tensor component suffix ---*/
+
+      fieldname.erase(fieldname.end()-3,fieldname.end());
+
+      AddDataArray(VTKDatatype::FLOAT32, fieldname, NTENSOR, myPoint*NTENSOR, GlobalPoint*NTENSOR);
+
+    } else if (output_variable && isVector) {
 
       /*--- Adjust the string name to remove the leading "X-" ---*/
 
@@ -252,33 +309,111 @@ void CParaviewXMLFileWriter::WriteData(string val_filename){
   VarCounter = varStart;
   for (iField = varStart; iField < fieldNames.size(); iField++) {
 
-    /*--- Check whether this field is a vector or scalar. ---*/
+    /*--- Check whether this field is a vector, tensor, or scalar. ---*/
+    /*--- Check tensor components first (longer patterns) to avoid double counting ---*/
 
-    bool output_variable = true, isVector = false;
-    size_t found = fieldNames[iField].find("_x");
+    bool output_variable = true, isVector = false, isTensor = false;
+
+    // Check for tensor components first (_xx, _xy, _yy, _xz, _yz, _zz)
+    size_t found = fieldNames[iField].find("_xx");
     if (found!=string::npos) {
       output_variable = true;
-      isVector        = true;
+      isTensor        = true;
+    }
+    found = fieldNames[iField].find("_xy");
+    if (found!=string::npos) {
+      /*--- We have found a tensor, so skip the XY component. ---*/
+      output_variable = false;
+      isTensor        = true;
+      VarCounter++;
     }
-    found = fieldNames[iField].find("_y");
+    found = fieldNames[iField].find("_yy");
     if (found!=string::npos) {
-      /*--- We have found a vector, so skip the Y component. ---*/
+      /*--- We have found a tensor, so skip the YY component. ---*/
       output_variable = false;
+      isTensor        = true;
       VarCounter++;
     }
-    found = fieldNames[iField].find("_z");
+    found = fieldNames[iField].find("_xz");
     if (found!=string::npos) {
-      /*--- We have found a vector, so skip the Z component. ---*/
+      /*--- We have found a tensor, so skip the XZ component. ---*/
       output_variable = false;
+      isTensor        = true;
+      VarCounter++;
+    }
+    found = fieldNames[iField].find("_yz");
+    if (found!=string::npos) {
+      /*--- We have found a tensor, so skip the YZ component. ---*/
+      output_variable = false;
+      isTensor        = true;
+      VarCounter++;
+    }
+    found = fieldNames[iField].find("_zz");
+    if (found!=string::npos) {
+      /*--- We have found a tensor, so skip the ZZ component. ---*/
+      output_variable = false;
+      isTensor        = true;
       VarCounter++;
     }
 
-    /*--- Write the point data as an <X,Y,Z> vector or a scalar. ---*/
+    // Check for vector components only if not a tensor (_x, _y, _z)
+    if (!isTensor) {
+      found = fieldNames[iField].find("_x");
+      if (found!=string::npos) {
+        output_variable = true;
+        isVector        = true;
+      }
+      found = fieldNames[iField].find("_y");
+      if (found!=string::npos) {
+        /*--- We have found a vector, so skip the Y component. ---*/
+        output_variable = false;
+        VarCounter++;
+      }
+      found = fieldNames[iField].find("_z");
+      if (found!=string::npos) {
+        /*--- We have found a vector, so skip the Z component. ---*/
+        output_variable = false;
+        VarCounter++;
+      }
+    }
+
+    /*--- Write the point data as a tensor, <X,Y,Z> vector, or scalar. ---*/
+
+    if (output_variable && isTensor) {
+
+      /*--- Resize buffer to accommodate tensor data ---*/
+      dataBufferFloat.resize(myPoint*NTENSOR);
+
+      for (iPoint = 0; iPoint < myPoint; iPoint++) {
+        dataBufferFloat[iPoint*NTENSOR + 0] = (float)dataSorter->GetData(VarCounter+0,iPoint); // XX
+        dataBufferFloat[iPoint*NTENSOR + 1] = (float)dataSorter->GetData(VarCounter+2,iPoint); // YY
+        dataBufferFloat[iPoint*NTENSOR + 3] = (float)dataSorter->GetData(VarCounter+1,iPoint); // XY
+        if (nDim == 2) {
+          /*--- For 2D tensors, set the off-diagonal z-components to 0, and zz-component to 1 ---*/
+
+          dataBufferFloat[iPoint*NTENSOR + 2] = 1.0; // ZZ = 1
+          dataBufferFloat[iPoint*NTENSOR + 4] = 0.0; // YZ = 0
+          dataBufferFloat[iPoint*NTENSOR + 5] = 0.0; // XZ = 0
+        } else {
+          /*--- For 3D tensors, get the z-components ---*/
+          dataBufferFloat[iPoint*NTENSOR + 2] = (float)dataSorter->GetData(VarCounter+5,iPoint); // ZZ
+          dataBufferFloat[iPoint*NTENSOR + 4] = (float)dataSorter->GetData(VarCounter+4,iPoint); // YZ
+          dataBufferFloat[iPoint*NTENSOR + 5] = (float)dataSorter->GetData(VarCounter+3,iPoint); // XZ
+        }
+      }
+
+      WriteDataArray(dataBufferFloat.data(), VTKDatatype::FLOAT32, myPoint*NTENSOR, GlobalPoint*NTENSOR,
+                     dataSorter->GetnPointCumulative(rank)*NTENSOR);
 
-    if (output_variable && isVector) {
+      VarCounter++;
+
+    } else if (output_variable && isVector) {
 
       /*--- Load up the buffer for writing this rank's vector data. ---*/
 
+      /*--- Resize buffer back to coordinate size ---*/
+      dataBufferFloat.resize(myPoint*NCOORDS);
+
       float val = 0.0;
       for (iPoint = 0; iPoint < myPoint; iPoint++) {
         for (iDim = 0; iDim < NCOORDS; iDim++) {
@@ -298,6 +433,8 @@ void CParaviewXMLFileWriter::WriteData(string val_filename){
 
     } else if (output_variable) {
 
+      /*--- For scalar data, ensure buffer is properly sized ---*/
+      dataBufferFloat.resize(myPoint);
 
       /*--- For now, create a temp 1D buffer to load up the data for writing.
        This will be replaced with a derived data type most likely. ---*/
diff --git a/SU2_CFD/src/solvers/CBaselineSolver.cpp b/SU2_CFD/src/solvers/CBaselineSolver.cpp
index fa89ee5e4830..ef6703f7efb2 100644
--- a/SU2_CFD/src/solvers/CBaselineSolver.cpp
+++ b/SU2_CFD/src/solvers/CBaselineSolver.cpp
@@ -40,6 +40,7 @@ CBaselineSolver::CBaselineSolver(CGeometry *geometry, CConfig *config) {
   /*--- Define geometry constants in the solver structure ---*/
 
   nDim = geometry->GetnDim();
+  nSymMat = 3 * (nDim - 1);
 
   /*--- Routines to access the number of variables and string names. ---*/
 
diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp
index 8bae5e6ca36f..2e8a1b3a0aaa 100644
--- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp
+++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp
@@ -41,6 +41,7 @@ CDiscAdjSolver::CDiscAdjSolver(CGeometry *geometry, CConfig *config, CSolver *di
 
   nVar = direct_solver->GetnVar();
   nDim = geometry->GetnDim();
+  nSymMat = 3 * (nDim - 1);
 
   nMarker      = config->GetnMarker_All();
   nPoint       = geometry->GetnPoint();
diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp
index dcca14edd19d..ebb02534b492 100644
--- a/SU2_CFD/src/solvers/CEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CEulerSolver.cpp
@@ -114,6 +114,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config,
    ---*/
 
   nDim = geometry->GetnDim();
+  nSymMat = 3 * (nDim - 1);
 
   nVar = nDim + 2;
   nPrimVar = nDim + 9;
diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
index 3df1ffbe3dd8..53d672829298 100644
--- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
@@ -104,6 +104,7 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned
    * Incompressible flow, primitive variables (P, vx, vy, vz, T, rho, beta, lamMu, EddyMu, Kt_eff, Cp, Cv) ---*/
 
   nDim = geometry->GetnDim();
+  nSymMat = 3 * (nDim - 1);
 
   /*--- Make sure to align the sizes with the constructor of CIncEulerVariable. ---*/
   nVar = nDim + 2;
diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp
index 65f3ea81f123..1497f4fe0f70 100644
--- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp
@@ -92,6 +92,7 @@ CNEMOEulerSolver::CNEMOEulerSolver(CGeometry *geometry, CConfig *config,
   nSpecies     = config->GetnSpecies();
   nMarker      = config->GetnMarker_All();
   nDim         = geometry->GetnDim();
+  nSymMat      = 3 * (nDim - 1);
   nPoint       = geometry->GetnPoint();
   nPointDomain = geometry->GetnPointDomain();
 
diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index a9c3875406bf..58fc19bd293f 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/SU2_CFD/src/solvers/CSolver.cpp
@@ -29,6 +29,7 @@
 #include "../../include/solvers/CSolver.hpp"
 #include "../../include/gradients/computeGradientsGreenGauss.hpp"
 #include "../../include/gradients/computeGradientsLeastSquares.hpp"
+#include "../../include/gradients/computeMetrics.hpp"
 #include "../../include/limiters/computeLimiters.hpp"
 #include "../../../Common/include/toolboxes/MMS/CIncTGVSolution.hpp"
 #include "../../../Common/include/toolboxes/MMS/CInviscidVortexSolution.hpp"
@@ -1379,6 +1380,14 @@ void CSolver::GetCommCountAndType(const CConfig* config,
       COUNT_PER_POINT  = nVar;
       MPI_TYPE         = COMM_TYPE_DOUBLE;
       break;
+    case MPI_QUANTITIES::GRADIENT_ADAPT:
+      COUNT_PER_POINT  = config->GetnMetric_Sensor()*nDim;
+      MPI_TYPE         = COMM_TYPE_DOUBLE;
+      break;
+    case MPI_QUANTITIES::HESSIAN:
+      COUNT_PER_POINT  = config->GetnMetric_Sensor()*nSymMat;
+      MPI_TYPE         = COMM_TYPE_DOUBLE;
+      break;
     default:
       SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.",
                      CURRENT_FUNCTION);
@@ -1393,6 +1402,8 @@ namespace CommHelpers {
       case MPI_QUANTITIES::PRIMITIVE_GRADIENT: return nodes->GetGradient_Primitive();
       case MPI_QUANTITIES::PRIMITIVE_GRAD_REC: return nodes->GetGradient_Reconstruction();
       case MPI_QUANTITIES::AUXVAR_GRADIENT: return nodes->GetAuxVarGradient();
+      case MPI_QUANTITIES::GRADIENT_ADAPT: return nodes->GetGradient_Adapt();
+      case MPI_QUANTITIES::HESSIAN: return nodes->GetHessian();
       default: return nodes->GetGradient();
     }
   }
@@ -1409,7 +1420,7 @@ void CSolver::InitiateComms(CGeometry *geometry,
 
   /*--- Local variables ---*/
 
-  unsigned short iVar, iDim;
+  unsigned short iVar, iDim, iMat;
   unsigned short COUNT_PER_POINT = 0;
   unsigned short MPI_TYPE        = 0;
 
@@ -1435,6 +1446,7 @@ void CSolver::InitiateComms(CGeometry *geometry,
   /*--- Handle the different types of gradient and limiter. ---*/
 
   const auto nVarGrad = COUNT_PER_POINT / nDim;
+  const auto nVarHess = COUNT_PER_POINT / nSymMat;
   auto& gradient = CommHelpers::selectGradient(base_nodes, commType);
   auto& limiter = CommHelpers::selectLimiter(base_nodes, commType);
 
@@ -1503,10 +1515,16 @@ void CSolver::InitiateComms(CGeometry *geometry,
           case MPI_QUANTITIES::SOLUTION_GRAD_REC:
           case MPI_QUANTITIES::PRIMITIVE_GRAD_REC:
           case MPI_QUANTITIES::AUXVAR_GRADIENT:
+          case MPI_QUANTITIES::GRADIENT_ADAPT:
             for (iVar = 0; iVar < nVarGrad; iVar++)
               for (iDim = 0; iDim < nDim; iDim++)
                 bufDSend[buf_offset+iVar*nDim+iDim] = gradient(iPoint, iVar, iDim);
             break;
+          case MPI_QUANTITIES::HESSIAN:
+            for (iVar = 0; iVar < nVarHess; iVar++)
+              for (iMat = 0; iMat < nSymMat; iMat++)
+                bufDSend[buf_offset+iVar*nSymMat+iMat] = gradient(iPoint, iVar, iMat);
+            break;
           case MPI_QUANTITIES::SOLUTION_FEA:
             for (iVar = 0; iVar < nVar; iVar++) {
               bufDSend[buf_offset+iVar] = base_nodes->GetSolution(iPoint, iVar);
@@ -1551,7 +1569,7 @@ void CSolver::CompleteComms(CGeometry *geometry,
 
   /*--- Local variables ---*/
 
-  unsigned short iDim, iVar;
+  unsigned short iDim, iVar, iMat;
   unsigned long iPoint, iRecv, nRecv, msg_offset, buf_offset;
   unsigned short COUNT_PER_POINT = 0;
   unsigned short MPI_TYPE = 0;
@@ -1572,6 +1590,7 @@ void CSolver::CompleteComms(CGeometry *geometry,
   /*--- Handle the different types of gradient and limiter. ---*/
 
   const auto nVarGrad = COUNT_PER_POINT / nDim;
+  const auto nVarHess = COUNT_PER_POINT / nSymMat;
   auto& gradient = CommHelpers::selectGradient(base_nodes, commType);
   auto& limiter = CommHelpers::selectLimiter(base_nodes, commType);
 
@@ -1651,10 +1670,16 @@ void CSolver::CompleteComms(CGeometry *geometry,
           case MPI_QUANTITIES::SOLUTION_GRAD_REC:
           case MPI_QUANTITIES::PRIMITIVE_GRAD_REC:
           case MPI_QUANTITIES::AUXVAR_GRADIENT:
+          case MPI_QUANTITIES::GRADIENT_ADAPT:
             for (iVar = 0; iVar < nVarGrad; iVar++)
               for (iDim = 0; iDim < nDim; iDim++)
                 gradient(iPoint,iVar,iDim) = bufDRecv[buf_offset+iVar*nDim+iDim];
             break;
+          case MPI_QUANTITIES::HESSIAN:
+            for (iVar = 0; iVar < nVarHess; iVar++)
+              for (iMat = 0; iMat < nSymMat; iMat++)
+                gradient(iPoint, iVar, iMat) = bufDRecv[buf_offset+iVar*nSymMat+iMat];
+            break;
           case MPI_QUANTITIES::SOLUTION_FEA:
             for (iVar = 0; iVar < nVar; iVar++) {
               base_nodes->SetSolution(iPoint, iVar, bufDRecv[buf_offset+iVar]);
@@ -2165,6 +2190,20 @@ void CSolver::SetSolution_Gradient_LS(CGeometry *geometry, const CConfig *config
   computeGradientsLeastSquares(this, comm, commPer, *geometry, *config, weighted, solution, 0, nVar, idxVel, gradient, rmatrix);
 }
 
+void CSolver::SetHessian_GG(CGeometry *geometry, const CConfig *config, short idxVel, const unsigned short Kind_Solver) {
+  const auto& solution = base_nodes->GetPrimitive_Adapt();
+  auto& gradient = base_nodes->GetGradient_Adapt();
+  auto nHess = config->GetnMetric_Sensor();
+
+  computeGradientsGreenGauss(this, MPI_QUANTITIES::GRADIENT_ADAPT, PERIODIC_GRAD_ADAPT,
+                             *geometry, *config, solution, 0, nHess, idxVel, gradient);
+
+  auto& hessian = base_nodes->GetHessian();
+
+  computeHessiansGreenGauss(this, MPI_QUANTITIES::HESSIAN, PERIODIC_HESSIAN,
+                            *geometry, *config, gradient, 0, nHess, idxVel, hessian);
+}
+
 void CSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *config) {
 
   /*--- Loop domain points. ---*/
@@ -4346,3 +4385,79 @@ void CSolver::SavelibROM(CGeometry *geometry, CConfig *config, bool converged) {
 #endif
 
 }
+
+
+void CSolver::ComputeMetric(CSolver **solver, CGeometry *geometry, const CConfig *config) {
+  /*--- TODO: - goal-oriented metric ---*/
+  /*---       - metric intersection  ---*/
+  const unsigned long nPointDomain = geometry->GetnPointDomain();
+  const unsigned short nSensor = config->GetnMetric_Sensor();
+
+  const bool normalize = (config->GetNormalize_Metric());
+
+
+  const unsigned long time_iter = config->GetTimeIter();
+  const bool steady = (config->GetTime_Marching() == TIME_MARCHING::STEADY);
+  const bool time_stepping = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) ||
+                             (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) ||
+                             (config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING);
+  const bool is_last_iter = (time_iter == config->GetnTime_Iter() - 1) || (steady);
+
+  /*--- Integrate and normalize the metric tensor field ---*/
+  vector<double> integrals;
+  for (auto iSensor = 0u; iSensor < nSensor; ++iSensor) {
+    SU2_OMP_MASTER
+    /*--- Make the Hessian eigenvalues positive definite, and add to the metric tensor ---*/
+    auto& hessians = base_nodes->GetHessian();
+    setPositiveDefiniteMetrics<su2double, tensor::hessian>(*geometry, *config, iSensor, hessians);
+    AddMetrics(solver, geometry, config, iSensor);
+
+    /*--- Integrate metric field on the last iteration (the end of the simulation if steady) ---*/
+    auto& metrics = base_nodes->GetMetric();
+    double integral = 0.0;
+    if (is_last_iter)
+      integral = integrateMetrics<double>(*geometry, *config, iSensor, metrics);
+
+    /*--- Normalize the metric field for steady simulations, or if requested for unsteady ---*/
+    if (steady || (normalize && is_last_iter))
+      normalizeMetrics<double, tensor::metric>(*geometry, *config, iSensor, integral, metrics);
+
+    /*--- Store the integral to be written ---*/
+    if (is_last_iter) {
+      integrals.push_back(integral);
+      if (rank == MASTER_NODE) {
+        cout << "Global metric normalization integral for sensor ";
+        cout << config->GetMetric_SensorString(iSensor) << ": " << integral << endl;
+      }
+    }
+    END_SU2_OMP_MASTER
+    SU2_OMP_BARRIER
+  }
+}
+
+void CSolver::AddMetrics(CSolver **solver, const CGeometry*geometry, const CConfig *config,
+                         const unsigned short iSensor) {
+  /*--- TODO: - goal-oriented metric ---*/
+  /*---       - metric intersection  ---*/
+  auto varFlo = solver[FLOW_SOL]->GetNodes();
+
+  const unsigned long nPointDomain = geometry->GetnPointDomain();
+  const unsigned short nSymMat = 3*(nDim-1);
+
+  const unsigned long time_iter = config->GetTimeIter();
+  const bool time_stepping = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) ||
+                             (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) ||
+                             (config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING);
+  const bool is_first_iter = (time_iter == 0);
+  const bool is_last_iter = (time_iter == config->GetnTime_Iter() - 1);
+
+  double coeff = (time_stepping && (is_first_iter || is_last_iter))? 0.5 : 1.0;
+  if (time_stepping) coeff *= SU2_TYPE::GetValue(config->GetTime_Step());
+
+  for(auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) {
+    for (auto iMat = 0; iMat < nSymMat; ++iMat) {
+      double hess = SU2_TYPE::GetValue(varFlo->GetHessian(iPoint, iSensor, iMat));
+      varFlo->AddMetric(iPoint, iMat, coeff * hess);
+    }
+  }
+}
diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp
index 7649165d0be3..0124d6e887b5 100644
--- a/SU2_CFD/src/solvers/CTurbSASolver.cpp
+++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp
@@ -53,6 +53,7 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor
   /*--- Define geometry constants in the solver structure ---*/
 
   nDim = geometry->GetnDim();
+  nSymMat = 3 * (nDim - 1);
 
   /*--- Single grid simulation ---*/
 
@@ -301,7 +302,7 @@ void CTurbSASolver::Viscous_Residual(const unsigned long iEdge, const CGeometry*
     /*--- Roughness heights. ---*/
     numerics->SetRoughness(geometry->nodes->GetRoughnessHeight(iPoint), geometry->nodes->GetRoughnessHeight(jPoint));
   };
-  
+
   /*--- Now instantiate the generic non-conservative implementation with the functor above. ---*/
   Viscous_Residual_NonCons(iEdge, geometry, solver_container, numerics, config, SolverSpecificNumerics);
 
diff --git a/SU2_CFD/src/variables/CFlowVariable.cpp b/SU2_CFD/src/variables/CFlowVariable.cpp
index 6d07949d79f2..2723543cee21 100644
--- a/SU2_CFD/src/variables/CFlowVariable.cpp
+++ b/SU2_CFD/src/variables/CFlowVariable.cpp
@@ -97,6 +97,13 @@ CFlowVariable::CFlowVariable(unsigned long npoint, unsigned long ndim, unsigned
   if (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE) {
     HB_Source.resize(nPoint, nVar) = su2double(0.0);
   }
+
+  if (config->GetCompute_Metric()) {
+    unsigned short nHess = config->GetnMetric_Sensor();
+    unsigned short nSymMat = 3 * (nDim - 1);
+    Primitive_Adapt.resize(nPoint, nHess) = su2double(0.0);
+    Metric.resize(nPoint, nSymMat) = 0.0;
+  }
 }
 
 void CFlowVariable::SetSolution_New() {
diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp
index 111a2a17775d..7c4ba2a70c16 100644
--- a/SU2_CFD/src/variables/CVariable.cpp
+++ b/SU2_CFD/src/variables/CVariable.cpp
@@ -51,6 +51,7 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva
   nPoint = npoint;
   nDim = ndim;
   nVar = nvar;
+  nSymMat = 3 * (nDim - 1);
 
   /*--- Allocate fields common to all problems. Do not allocate fields
    that are specific to one solver, i.e. not common, in this class. ---*/
@@ -79,6 +80,13 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva
 
   if (config->GetMultizone_Problem())
     Solution_BGS_k.resize(nPoint,nVar) = su2double(0.0);
+
+  /*--- Gradient and Hessian for anisotropic metric ---*/
+  if (config->GetCompute_Metric()) {
+    unsigned short nHess = config->GetnMetric_Sensor();
+    Gradient_Adapt.resize(nPoint,nHess,nDim,0.0);
+    Hessian.resize(nPoint,nHess,nSymMat,0.0);
+  }
 }
 
 void CVariable::Set_OldSolution() {