diff --git a/src/optimisation/adjointOptimisation/adjoint/Make/options b/src/optimisation/adjointOptimisation/adjoint/Make/options
index 44dec84580..f3041856b6 100644
--- a/src/optimisation/adjointOptimisation/adjoint/Make/options
+++ b/src/optimisation/adjointOptimisation/adjoint/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
@@ -14,6 +15,7 @@ EXE_INC = \
 
 LIB_LIBS = \
     -lfiniteVolume \
+    -lfiniteArea \
     -lmeshTools \
     -lsurfMesh \
     -lsampling \
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C
index 90bcce43aa..e17c5365a8 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2020 PCOpt/NTUA
-    Copyright (C) 2013-2020 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,11 @@ License
 #include "PrimitivePatchInterpolation.H"
 #include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
+#include "faCFD.H"
+#include "fixedValueFaPatchFieldsFwd.H"
+#include "fixedValueFaPatchFields.H"
+#include "zeroGradientFaPatchFieldsFwd.H"
+#include "zeroGradientFaPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -42,7 +47,7 @@ namespace incompressible
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-defineTypeNameAndDebug(sensitivitySurface, 0);
+defineTypeNameAndDebug(sensitivitySurface, 1);
 addToRunTimeSelectionTable
 (
     adjointSensitivity,
@@ -214,6 +219,188 @@ void sensitivitySurface::setSuffixName()
 }
 
 
+void sensitivitySurface::smoothSensitivities()
+{
+    // Read in parameters
+    const label iters(dict().getOrDefault<label>("iters", 500));
+    const scalar tolerance(dict().getOrDefault<scalar>("tolerance", 1.e-06));
+    autoPtr<faMesh> aMeshPtr(nullptr);
+
+    IOobject faceLabels
+    (
+        "faceLabels",
+        mesh_.time().findInstance
+        (
+            mesh_.dbDir()/faMesh::meshSubDir, 
+            "faceLabels", 
+            IOobject::READ_IF_PRESENT
+        ),
+        faMesh::meshSubDir,
+        mesh_,
+        IOobject::READ_IF_PRESENT,
+        IOobject::NO_WRITE
+    );
+
+    // If the faMesh already exists, read it
+    if (faceLabels.typeHeaderOk<labelIOList>(false))
+    {
+        Info<< "Reading the already constructed faMesh" << endl;
+        aMeshPtr.reset(new faMesh(mesh_));
+    }
+    else
+    {
+        // Dictionary used to construct the faMesh
+        dictionary faMeshDefinition;
+
+        IOobject faMeshDefinitionDict
+        (
+            "faMeshDefinition",
+            mesh_.time().caseSystem(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+        
+        // If the faMeshDefinitionDict exists, use it to construct the mesh
+        if (faMeshDefinitionDict.typeHeaderOk<IOdictionary>(false))
+        {
+            Info<< "Reading faMeshDefinition from system " << endl;
+            faMeshDefinition = IOdictionary(faMeshDefinitionDict);
+        }
+        // Otherwise, faMesh is generated from all patches on which we compute
+        // sensitivities
+        else
+        {
+            Info<< "Constructing faMeshDefinition from sensitivity patches" 
+                << endl;
+            wordList polyMeshPatches(sensitivityPatchIDs_.size());
+            label i(0);
+            for (const label patchID : sensitivityPatchIDs_)
+            {
+                polyMeshPatches[i++] = mesh_.boundary()[patchID].name();    
+            }
+            faMeshDefinition.add<wordList>("polyMeshPatches", polyMeshPatches);
+            dictionary& boundary = faMeshDefinition.subDictOrAdd("boundary");
+        }
+
+        // Construct faMesh
+        aMeshPtr.reset(new faMesh(mesh_, faMeshDefinition));
+    }
+    faMesh& aMesh = aMeshPtr.ref();
+
+    // Physical radius of the smoothing, provided either directly or computed
+    // based on the average 'length' of boundary faces
+    const scalar Rphysical
+        (dict().getOrDefault<scalar>("radius", computeRadius(aMesh)));
+    DebugInfo
+        << "Physical radius of the sensitivity smoothing " 
+        << Rphysical << nl << endl;
+
+    // Radius used as the diffusivity in the Helmholtz filter, computed as a 
+    // function of the physical radius
+    const dimensionedScalar RpdeSqr
+    (
+        "RpdeSqr", dimArea, sqr(Rphysical/(2.*::sqrt(3.)))
+    );
+
+    dimensionedScalar one("1", dimless, 1.);
+
+    // Mapping engine
+    volSurfaceMapping vsm(aMesh);
+
+    // Source term in faMatrix needs to be an areaField
+    areaScalarField sens
+    (
+        IOobject
+        (
+            "sens",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        aMesh,
+        dimensionedScalar(dimless, Zero),
+        zeroGradientFaPatchField<scalar>::typeName
+    );
+
+    // Copy sensitivities to area field
+    sens.primitiveFieldRef() =
+        vsm.mapToSurface<scalar>(wallFaceSensNormalPtr_());
+
+    // Initialisation of the smoothed sensitivities field based on the original
+    // sensitivities
+    areaScalarField smoothedSens("smoothedSens", sens);
+    for (label iter = 0; iter < iters; ++iter)
+    {
+        Info<< "Sensitivity smoothing iteration " << iter << endl;
+
+        faScalarMatrix smoothEqn
+        (
+            fam::Sp(one, smoothedSens)
+          - fam::laplacian(RpdeSqr, smoothedSens)
+         ==
+            sens
+        );
+
+        smoothEqn.relax();
+
+        const scalar residual(mag(smoothEqn.solve().initialResidual()));
+
+        DebugInfo
+            << "Max smoothSens " << gMax(mag(smoothedSens)()) << endl;
+
+        // Print execution time
+        mesh_.time().printExecutionTime(Info);
+
+        // Check convergence
+        if (residual < tolerance)
+        {
+            Info<< "\n***Reached smoothing equation convergence limit, "
+                   "iteration " << iter << "***\n\n";
+            break;
+        }
+    }
+
+    // Field used to write the smoothed sensitivity field to file
+    volScalarField volSmoothedSens
+    (
+        IOobject
+        (
+            "smoothedSurfaceSens" + surfaceFieldSuffix_,
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar(dimless, Zero)
+    );
+
+    // Transfer result back to volField and write
+    vsm.mapToVolume(smoothedSens, volSmoothedSens.boundaryFieldRef());
+    volSmoothedSens.write();
+}
+
+
+scalar sensitivitySurface::computeRadius(const faMesh& aMesh)
+{
+    scalar averageArea(gAverage(aMesh.S().field()));
+    const Vector<label>& geometricD = mesh_.geometricD();
+    const boundBox& bounds = mesh_.bounds();
+    forAll(geometricD, iDir)
+    {
+        if (geometricD[iDir] == -1)
+        {
+            averageArea /= bounds.span()[iDir];
+        }
+    }
+    scalar mult = dict().getOrDefault<scalar>("meanRadiusMultiplier", 10);
+
+    return mult*pow(averageArea, scalar(1)/scalar(mesh_.nGeometricD() - 1));
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 sensitivitySurface::sensitivitySurface
@@ -244,6 +431,7 @@ sensitivitySurface::sensitivitySurface
     includeMeshMovement_(false),
     includeObjective_(false),
     writeGeometricInfo_(false),
+    smoothSensitivities_(false),
     eikonalSolver_(nullptr),
     meshMovementSolver_(nullptr),
 
@@ -346,6 +534,8 @@ void sensitivitySurface::read()
         dict().getOrDefault<bool>("includeObjectiveContribution", true);
     writeGeometricInfo_ =
         dict().getOrDefault<bool>("writeGeometricInfo", false);
+    smoothSensitivities_ =
+        dict().getOrDefault<bool>("smoothSensitivities", false);
 
     // Allocate new solvers if necessary
     if (includeDistance_ && !eikonalSolver_)
@@ -700,6 +890,12 @@ void sensitivitySurface::assembleSensitivities()
         }
         nPassedFaces += patch.size();
     }
+
+    // Smooth sensitivities if needed
+    if (smoothSensitivities_)
+    {
+        smoothSensitivities();
+    }
 }
 
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.H
index ec070cb0da..baff1c1ecf 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.H
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.H
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2020 PCOpt/NTUA
-    Copyright (C) 2013-2020 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    Copyright (C) 2013-2021 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,27 @@ Class
 Description
     Calculation of adjoint based sensitivities at wall faces
 
+    Reference:
+    \verbatim
+        The computation of the sensitivity derivatives  follows the (E-)SI
+        formulation of
+
+            Kavvadias, I. S., Papoutsis-Kiachagias, E. M., 
+            Giannakoglou, K. C. (2015).
+            On the proper treatment of grid sensitivities in continuous adjoint
+            methods for shape optimization.
+            Journal of computational physics, 301, 1-18.
+            https://doi.org/10.1016/j.jcp.2015.08.012
+
+        whereas their smoothing based on the computation of the 'Sobolev
+        gradient' is derived from
+    
+            Vassberg J. C., Jameson A. (2006).
+            Aerodynamic Shape Optimization Part I: Theoretical Background.
+            VKI Lecture Series, Introduction to Optimization and
+            Multidisciplinary Design, Brussels, Belgium, 8 March, 2006.
+    \endverbatim
+
 SourceFiles
     sensitivitySurface.C
 
@@ -44,6 +65,7 @@ SourceFiles
 #include "adjointEikonalSolverIncompressible.H"
 #include "adjointMeshMovementSolverIncompressible.H"
 #include "deltaBoundary.H"
+#include "faMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -97,6 +119,10 @@ protected:
         //- Write geometric info for use by external programs
         bool writeGeometricInfo_;
 
+        //- Smooth sensitivity derivatives based on the computation of the 
+        //- 'Sobolev gradient'
+        bool smoothSensitivities_;
+
         autoPtr<adjointEikonalSolver> eikonalSolver_;
 
         autoPtr<adjointMeshMovementSolver> meshMovementSolver_;
@@ -116,6 +142,14 @@ protected:
         //- Set suffix name for sensitivity fields
         void setSuffixName();
 
+        //- Smooth sensitivity derivatives based on the computation of the 
+        //- 'Sobolev gradient'
+        void smoothSensitivities();
+
+        //- Compute the physical smoothing radius based on the average boundary
+        //- face 'length'
+        scalar computeRadius(const faMesh& aMesh);
+
 
 private:
 
@@ -177,7 +211,7 @@ public:
        //- Write sensitivity maps
        virtual void write(const word& baseName = word::null);
 
-       // Inline geters and setters
+       // Inline getters and setters
 
            //- Get access to the includeObjective bool
            inline bool getIncludeObjective() const;