diff --git a/src/optimisation/adjointOptimisation/adjoint/Make/files b/src/optimisation/adjointOptimisation/adjoint/Make/files
index d79ad9af80ae8c317d0c467303f52d8f8f2073f9..eb253511fbf7f61a466618f75973926bfe0762cc 100644
--- a/src/optimisation/adjointOptimisation/adjoint/Make/files
+++ b/src/optimisation/adjointOptimisation/adjoint/Make/files
@@ -126,6 +126,7 @@ optimisation/adjointSensitivity/sensitivity/sensitivity.C
 optimisation/adjointSensitivity/shapeSensitivitiesBase/shapeSensitivitiesBase.C
 incoSens=optimisation/adjointSensitivity/incompressible
 $(incoSens)/adjointSensitivity/adjointSensitivityIncompressible.C
+$(incoSens)/shapeSensitivities/shapeSensitivitiesIncompressible.C
 $(incoSens)/adjointEikonalSolver/adjointEikonalSolverIncompressible.C
 $(incoSens)/adjointMeshMovementSolver/adjointMeshMovementSolverIncompressible.C
 $(incoSens)/sensitivitySurface/sensitivitySurfaceIncompressible.C
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/FIBase/FIBaseIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/FIBase/FIBaseIncompressible.C
index ca4fc90810ac7253f50918b550b4bd214f2d5841..4cc2993e2aff48e28c3d24f3e3168949b35ff5f3 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/FIBase/FIBaseIncompressible.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/FIBase/FIBaseIncompressible.C
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -84,7 +84,7 @@ FIBase::FIBase
     fv::optionAdjointList& fvOptionsAdjoint
 )
 :
-    adjointSensitivity
+    shapeSensitivities
     (
         mesh,
         dict,
@@ -93,7 +93,6 @@ FIBase::FIBase
         objectiveManager,
         fvOptionsAdjoint
     ),
-    shapeSensitivitiesBase(mesh, dict),
     gradDxDbMult_
     (
         IOobject
@@ -109,9 +108,6 @@ FIBase::FIBase
     ),
     divDxDbMult_(mesh_.nCells(), Zero),
     optionsDxDbMult_(mesh_.nCells(), Zero),
-    dSfdbMult_(createZeroBoundaryPtr<vector>(mesh_)),
-    dnfdbMult_(createZeroBoundaryPtr<vector>(mesh_)),
-    dxdbDirectMult_(createZeroBoundaryPtr<vector>(mesh_)),
 
     includeDistance_(false),
     eikonalSolver_(nullptr)
@@ -165,18 +161,11 @@ void FIBase::accumulateIntegrand(const scalar dt)
     }
 
     // Accumulate direct sensitivities
-    for (const label patchI : sensitivityPatchIDs_)
-    {
-        const scalarField magSfDt(mesh_.boundary()[patchI].magSf()*dt);
-        for (objective& func : functions)
-        {
-            const scalar wei = func.weight();
-            dSfdbMult_()[patchI] += wei*func.dSdbMultiplier(patchI)*dt;
-            dnfdbMult_()[patchI] += wei*func.dndbMultiplier(patchI)*magSfDt;
-            dxdbDirectMult_()[patchI] +=
-                wei*func.dxdbDirectMultiplier(patchI)*magSfDt;
-        }
-    }
+    accumulateDirectSensitivityIntegrand(dt);
+
+    // Accumulate sensitivities due to boundary conditions
+    accumulateBCSensitivityIntegrand(dt);
+
 }
 
 
@@ -185,24 +174,13 @@ void FIBase::clearSensitivities()
     gradDxDbMult_ = dimensionedTensor(gradDxDbMult_.dimensions(), Zero);
     divDxDbMult_ = Zero;
     optionsDxDbMult_ = vector::zero;
-    dSfdbMult_() = vector::zero;
-    dnfdbMult_() = vector::zero;
-    dxdbDirectMult_() = vector::zero;
 
     if (includeDistance_)
     {
         eikonalSolver_->reset();
     }
 
-    adjointSensitivity::clearSensitivities();
-    shapeSensitivitiesBase::clear();
-}
-
-
-void FIBase::write(const word& baseName)
-{
-    adjointSensitivity::write(baseName);
-    shapeSensitivitiesBase::write();
+    shapeSensitivities::clearSensitivities();
 }
 
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/FIBase/FIBaseIncompressible.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/FIBase/FIBaseIncompressible.H
index e7f1120caeafff3e4f486c873633dd0284cead09..576493b72a43c3b63c13a291a7a618606572a8b3 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/FIBase/FIBaseIncompressible.H
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/FIBase/FIBaseIncompressible.H
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -40,7 +40,7 @@ SourceFiles
 #define FIBaseIncompressible_H
 
 #include "adjointSensitivityIncompressible.H"
-#include "shapeSensitivitiesBase.H"
+#include "shapeSensitivitiesIncompressible.H"
 #include "adjointEikonalSolverIncompressible.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -52,13 +52,12 @@ namespace incompressible
 {
 
 /*---------------------------------------------------------------------------*\
-                  Class FIBase Declaration
+                            Class FIBase Declaration
 \*---------------------------------------------------------------------------*/
 
 class FIBase
 :
-    public adjointSensitivity,
-    public shapeSensitivitiesBase
+    public shapeSensitivities
 {
 protected:
 
@@ -73,11 +72,6 @@ protected:
         //- dx/db multiplier coming from fvOptions
         vectorField optionsDxDbMult_;
 
-        //- Fields related to direct sensitivities
-        autoPtr<boundaryVectorField> dSfdbMult_;
-        autoPtr<boundaryVectorField> dnfdbMult_;
-        autoPtr<boundaryVectorField> dxdbDirectMult_;
-
         //- Include distance variation in sens computation
         bool includeDistance_;
 
@@ -139,9 +133,6 @@ public:
 
         //- Zero sensitivity fields and their constituents
         virtual void clearSensitivities();
-
-        //- Write sensitivity fields
-        virtual void write(const word& baseName = word::null);
 };
 
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/SIBase/SIBaseIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/SIBase/SIBaseIncompressible.C
index ed34fc5ccd78c0f47e8e2679c3233853d74fd36c..96465d521afff9a2bbd79587cac5696a05bd76c4 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/SIBase/SIBaseIncompressible.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/SIBase/SIBaseIncompressible.C
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -83,7 +83,7 @@ SIBase::SIBase
     fv::optionAdjointList& fvOptionsAdjoint
 )
 :
-    adjointSensitivity
+    shapeSensitivities
     (
         mesh,
         dict,
@@ -92,12 +92,11 @@ SIBase::SIBase
         objectiveManager,
         fvOptionsAdjoint
     ),
-    shapeSensitivitiesBase(mesh, dict),
     surfaceSensitivity_
     (
         mesh,
         // Ideally, subOrEmptyDict would be used.
-        // Since we need a recursive search in shapeSensitivitiesBase though
+        // Since we need a recursive search in shapeSensitivities though
         // and the dict returned by subOrEmptyDict (if found)
         // does not know its parent, optionalSubDict is used
         dict.optionalSubDict("surfaceSensitivities"),
@@ -106,9 +105,6 @@ SIBase::SIBase
         objectiveManager,
         fvOptionsAdjoint
     ),
-    dSfdbMult_(createZeroBoundaryPtr<vector>(mesh_)),
-    dnfdbMult_(createZeroBoundaryPtr<vector>(mesh_)),
-    dxdbDirectMult_(createZeroBoundaryPtr<vector>(mesh_)),
     includeObjective_(true)
 {
     read();
@@ -141,42 +137,18 @@ void SIBase::accumulateIntegrand(const scalar dt)
     // Accumulate direct sensitivities
     if (includeObjective_)
     {
-        PtrList<objective>& functions
-        (
-            objectiveManager_.getObjectiveFunctions()
-        );
-        for (const label patchI : sensitivityPatchIDs_)
-        {
-            const scalarField magSfDt(mesh_.boundary()[patchI].magSf()*dt);
-            for (objective& func : functions)
-            {
-                const scalar wei = func.weight();
-                dSfdbMult_()[patchI] += wei*func.dSdbMultiplier(patchI)*dt;
-                dnfdbMult_()[patchI] += wei*func.dndbMultiplier(patchI)*magSfDt;
-                dxdbDirectMult_()[patchI] +=
-                    wei*func.dxdbDirectMultiplier(patchI)*magSfDt;
-            }
-        }
+        accumulateDirectSensitivityIntegrand(dt);
     }
+
+    // Accumulate sensitivities due to boundary conditions
+    accumulateBCSensitivityIntegrand(dt);
 }
 
 
 void SIBase::clearSensitivities()
 {
     surfaceSensitivity_.clearSensitivities();
-    dSfdbMult_() = vector::zero;
-    dnfdbMult_() = vector::zero;
-    dxdbDirectMult_() = vector::zero;
-
-    adjointSensitivity::clearSensitivities();
-    shapeSensitivitiesBase::clear();
-}
-
-
-void SIBase::write(const word& baseName)
-{
-    adjointSensitivity::write(baseName);
-    shapeSensitivitiesBase::write();
+    shapeSensitivities::clearSensitivities();
 }
 
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/SIBase/SIBaseIncompressible.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/SIBase/SIBaseIncompressible.H
index 25de490626adfe3ed02bdb0f65c8c42d94c37a1b..ee97adac0556784848a0e37714391180c0649c3a 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/SIBase/SIBaseIncompressible.H
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/SIBase/SIBaseIncompressible.H
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -40,8 +40,7 @@ SourceFiles
 #ifndef SIBaseIncompressible_H
 #define SIBaseIncompressible_H
 
-#include "adjointSensitivityIncompressible.H"
-#include "shapeSensitivitiesBase.H"
+#include "shapeSensitivitiesIncompressible.H"
 #include "sensitivitySurfaceIncompressible.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -53,13 +52,12 @@ namespace incompressible
 {
 
 /*---------------------------------------------------------------------------*\
-                  Class SIBase Declaration
+                            Class SIBase Declaration
 \*---------------------------------------------------------------------------*/
 
 class SIBase
 :
-    public adjointSensitivity,
-    public shapeSensitivitiesBase
+    public shapeSensitivities
 {
 protected:
 
@@ -68,11 +66,10 @@ protected:
         //- Surface sensitivities
         sensitivitySurface surfaceSensitivity_;
 
-        //- Fields related to direct sensitivities
-        autoPtr<boundaryVectorField> dSfdbMult_;
-        autoPtr<boundaryVectorField> dnfdbMult_;
-        autoPtr<boundaryVectorField> dxdbDirectMult_;
-
+        //- Whether to include direct sensitivities or not.
+        //  Used to avoid double contributions from both here and the
+        //  sensitivitySurface object which might have already accounted for
+        //  them
         bool includeObjective_;
 
 
@@ -130,9 +127,6 @@ public:
 
         //- Zero sensitivity fields and their constituents
         virtual void clearSensitivities();
-
-        //- Write sensitivity fields.
-        virtual void write(const word& baseName = word::null);
 };
 
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezier/sensitivityBezierIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezier/sensitivityBezierIncompressible.C
index 4e1b1cfd5057c46dbf36b07e12d1bdce02e8a936..0f3235759245d5e98a961c24cf1add717ae48ec2 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezier/sensitivityBezierIncompressible.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezier/sensitivityBezierIncompressible.C
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -77,6 +77,7 @@ sensitivityBezier::sensitivityBezier
     dSdbSens_(Bezier_.nBezier(), Zero),
     dndbSens_(Bezier_.nBezier(), Zero),
     dxdbDirectSens_(Bezier_.nBezier(), Zero),
+    bcSens_(Bezier_.nBezier(), Zero),
     derivativesFolder_("optimisation"/type() + "Derivatives")
 {
     derivatives_ = scalarField(3*Bezier_.nBezier(), Zero);
@@ -137,9 +138,12 @@ void sensitivityBezier::assembleSensitivities()
                 dxdbDirectSens_[iCP] +=
                     gSum((dxdbDirectMult_()[patchI] & dxidXj));
             }
+
+            // Sensitivities from boundary conditions
+            bcSens_[iCP] += gSum(bcDxDbMult_()[patchI] & dxidXj);
         }
     }
-    sens_ = flowSens_ + dSdbSens_ + dndbSens_ + dxdbDirectSens_;
+    sens_ = flowSens_ + dSdbSens_ + dndbSens_ + dxdbDirectSens_ + bcSens_;
 
     // Transform sensitivities to scalarField in order to cooperate with
     // updateMethod
@@ -159,6 +163,7 @@ void sensitivityBezier::assembleSensitivities()
             dSdbSens_[cpI].x() = Zero;
             dndbSens_[cpI].x() = Zero;
             dxdbDirectSens_[cpI].x() = Zero;
+            bcSens_[cpI].x() = Zero;
         }
         if (confineYmovement[cpI])
         {
@@ -167,6 +172,7 @@ void sensitivityBezier::assembleSensitivities()
             dSdbSens_[cpI].y() = Zero;
             dndbSens_[cpI].y() = Zero;
             dxdbDirectSens_[cpI].y() = Zero;
+            bcSens_[cpI].y() = Zero;
         }
         if (confineZmovement[cpI])
         {
@@ -175,6 +181,7 @@ void sensitivityBezier::assembleSensitivities()
             dSdbSens_[cpI].z() = Zero;
             dndbSens_[cpI].z() = Zero;
             dxdbDirectSens_[cpI].z() = Zero;
+            bcSens_[cpI].z() = Zero;
         }
     }
 }
@@ -187,6 +194,7 @@ void sensitivityBezier::clearSensitivities()
     dSdbSens_ = Zero;
     dndbSens_ = Zero;
     dxdbDirectSens_ = Zero;
+    bcSens_ = Zero;
 
     SIBase::clearSensitivities();
 }
@@ -210,7 +218,8 @@ void sensitivityBezier::write(const word& baseName)
             << setw(width) << "flow" << " "
             << setw(width) << "dSdb" << " "
             << setw(width) << "dndb" << " "
-            << setw(width) << "dxdbDirect" << endl;
+            << setw(width) << "dxdbDirect" << " "
+            << setw(width) << "dvdb" << endl;
         label nDV = derivatives_.size();
         label nBezier = Bezier_.nBezier();
         const boolListList& confineMovement = Bezier_.confineMovement();
@@ -229,7 +238,8 @@ void sensitivityBezier::write(const word& baseName)
                     << setw(width) << flowSens_[iCP].component(idir) << " "
                     << setw(width) << dSdbSens_[iCP].component(idir) << " "
                     << setw(width) << dndbSens_[iCP].component(idir) << " "
-                    << setw(width) << dxdbDirectSens_[iCP].component(idir)
+                    << setw(width) << dxdbDirectSens_[iCP].component(idir) << " "
+                    << setw(width) << bcSens_[iCP].component(idir)
                     << endl;
             }
         }
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezier/sensitivityBezierIncompressible.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezier/sensitivityBezierIncompressible.H
index cd450e0e115df15c6491ab29366398662e37c5eb..ad5a99222b8d3fa3c9abe3b8d4a0a22d447f1531 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezier/sensitivityBezierIncompressible.H
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezier/sensitivityBezierIncompressible.H
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -78,6 +78,7 @@ protected:
         vectorField dSdbSens_;
         vectorField dndbSens_;
         vectorField dxdbDirectSens_;
+        vectorField bcSens_;
 
         fileName derivativesFolder_;
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezierFI/sensitivityBezierFIIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezierFI/sensitivityBezierFIIncompressible.C
index 72b6f425869491d35713562be89ed5c428fa607b..ca17f29794670a2d1775fe0988aff9a5e1fd6731 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezierFI/sensitivityBezierFIIncompressible.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezierFI/sensitivityBezierFIIncompressible.C
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -146,6 +146,7 @@ sensitivityBezierFI::sensitivityBezierFI
     dVdbSens_(3*Bezier_.nBezier(), Zero),
     distanceSens_(3*Bezier_.nBezier(), Zero),
     optionsSens_(3*Bezier_.nBezier(), Zero),
+    bcSens_(3*Bezier_.nBezier(), Zero),
 
     derivativesFolder_("optimisation"/type() + "Derivatives"),
 
@@ -249,6 +250,9 @@ void sensitivityBezierFI::assembleSensitivities()
                 const vectorField& dxdbFace = tdxdbFace();
                 dxdbDirectSens_[iDV] +=
                     gSum(dxdbDirectMult_()[patchI] & dxdbFace);
+
+                // Contribution from boundary conditions
+                bcSens_[iDV] += gSum(bcDxDbMult_()[patchI] & dxdbFace);
             }
 
             // Contribution from delta (V) / delta b
@@ -290,7 +294,8 @@ void sensitivityBezierFI::assembleSensitivities()
           + dxdbDirectSens_
           + dVdbSens_
           + distanceSens_
-          + optionsSens_;
+          + optionsSens_
+          + bcSens_;
     }
 }
 
@@ -304,6 +309,7 @@ void sensitivityBezierFI::clearSensitivities()
     dVdbSens_ = Zero;
     distanceSens_ = Zero;
     optionsSens_ = Zero;
+    bcSens_ = Zero;
 
     FIBase::clearSensitivities();
 }
@@ -329,7 +335,9 @@ void sensitivityBezierFI::write(const word& baseName)
             << setw(width)   << "dndb"       << " "
             << setw(width)   << "dxdbDirect" << " "
             << setw(width)   << "dVdb"       << " "
-            << setw(width)   << "distance"   << endl;
+            << setw(width)   << "distance"   << " "
+            << setw(width)   << "options"    << " "
+            << setw(width)   << "dvdb"       << endl;
         const label nDVs = derivatives_.size();
         const label nBezier = Bezier_.nBezier();
         const boolListList& confineMovement = Bezier_.confineMovement();
@@ -337,11 +345,14 @@ void sensitivityBezierFI::write(const word& baseName)
 
         for (label iDV = 0; iDV < nDVs; iDV++)
         {
-            label iCP = iDV%nBezier;
-            label idir = iDV/nBezier;
+            const label iCP(iDV%nBezier);
+            const label idir(iDV/nBezier);
             if (!confineMovement[idir][iCP])
             {
-                if (iDV!=lastActive + 1) derivFile << "\n";
+                if (iDV!=lastActive + 1)
+                {
+                    derivFile << "\n";
+                }
                 lastActive = iDV;
                 derivFile
                    << setw(widthDV) << iDV << " "
@@ -351,7 +362,9 @@ void sensitivityBezierFI::write(const word& baseName)
                    << setw(width) << dndbSens_[iDV] << " "
                    << setw(width) << dxdbDirectSens_[iDV] << " "
                    << setw(width) << dVdbSens_[iDV] << " "
-                   << setw(width) << distanceSens_[iDV] << endl;
+                   << setw(width) << distanceSens_[iDV] << " "
+                   << setw(width) << optionsSens_[iDV] << " "
+                   << setw(width) << bcSens_[iDV] << endl;
             }
         }
     }
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezierFI/sensitivityBezierFIIncompressible.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezierFI/sensitivityBezierFIIncompressible.H
index 382316b432bd29c7cedb03cc348b5f083598cabb..2ea45cb5fe3ad28fcfed46e19a4b2f8ef930634d 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezierFI/sensitivityBezierFIIncompressible.H
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityBezierFI/sensitivityBezierFIIncompressible.H
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -93,6 +93,9 @@ protected:
         //- Term depending on fvOptions
         scalarField optionsSens_;
 
+        //- Term depending on the differenation of boundary conditions
+        scalarField bcSens_;
+
         fileName derivativesFolder_;
 
         label meshMovementIters_;
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 3869c7dcf4605a0a4356153ff8b6261ce9d318f3..364c522f1570c1ee0c6712a1ef62b6c704999f48 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-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -707,7 +707,7 @@ void sensitivitySurface::clearSensitivities()
     }
     // Reset sensitivity fields
     adjointSensitivity::clearSensitivities();
-    shapeSensitivitiesBase::clear();
+    shapeSensitivitiesBase::clearSensitivities();
 }
 
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C
index f856d24072755853f1a16f1de546e4be2aacde6e..4716b65cb172ee8d5ba7d3683bc9dfdf886c6ed5 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -743,7 +743,7 @@ void sensitivitySurfacePoints::clearSensitivities()
 
     // Reset sensitivity fields
     adjointSensitivity::clearSensitivities();
-    shapeSensitivitiesBase::clear();
+    shapeSensitivitiesBase::clearSensitivities();
 }
 
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplines/sensitivityVolBSplinesIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplines/sensitivityVolBSplinesIncompressible.C
index 7e9e54a6352edccffad789022fe886a431fbed2f..66fb1553638fc33c40ab183bd49e3dd94a3314db 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplines/sensitivityVolBSplinesIncompressible.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplines/sensitivityVolBSplinesIncompressible.C
@@ -113,6 +113,30 @@ void sensitivityVolBSplines::computeObjectiveContributions()
 }
 
 
+void sensitivityVolBSplines::computeBCContributions()
+{
+    label passedCPs = 0;
+    PtrList<NURBS3DVolume>& boxes = volBSplinesBase_.boxesRef();
+    forAll(boxes, iNURB)
+    {
+        vectorField sensBcsDxDb =
+            boxes[iNURB].computeControlPointSensitivities
+            (
+                bcDxDbMult_(),
+                sensitivityPatchIDs_.toc()
+            );
+
+        // Transfer to global list
+        forAll(sensBcsDxDb, cpI)
+        {
+            bcSens_[passedCPs + cpI] = sensBcsDxDb[cpI];
+        }
+        passedCPs += sensBcsDxDb.size();
+    }
+    volBSplinesBase_.boundControlPointMovement(bcSens_);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 sensitivityVolBSplines::sensitivityVolBSplines
@@ -143,16 +167,18 @@ sensitivityVolBSplines::sensitivityVolBSplines
     dSdbSens_(0),
     dndbSens_(0),
     dxdbDirectSens_(0),
+    bcSens_(0),
 
     derivativesFolder_("optimisation"/type() + "Derivatives")
 {
     // No boundary field pointers need to be allocated
-    label nCPs = volBSplinesBase_.getTotalControlPointsNumber();
+    const label nCPs(volBSplinesBase_.getTotalControlPointsNumber());
     derivatives_ = scalarField(3*nCPs, Zero);
     flowSens_ = vectorField(nCPs, Zero);
     dSdbSens_ = vectorField(nCPs, Zero);
     dndbSens_ = vectorField(nCPs, Zero);
     dxdbDirectSens_ = vectorField(nCPs, Zero);
+    bcSens_ = vectorField(nCPs, Zero);
 
     // Create folder to store sensitivities
     mkDir(derivativesFolder_);
@@ -199,6 +225,8 @@ void sensitivityVolBSplines::assembleSensitivities()
     // false for sensitivityVolBSplines
     computeObjectiveContributions();
 
+    computeBCContributions();
+
     // Transform sensitivites to scalarField in order to cooperate with
     // updateMethod
     forAll(flowSens_, cpI)
@@ -207,17 +235,20 @@ void sensitivityVolBSplines::assembleSensitivities()
             flowSens_[cpI].x()
           + dSdbSens_[cpI].x()
           + dndbSens_[cpI].x()
-          + dxdbDirectSens_[cpI].x();
+          + dxdbDirectSens_[cpI].x()
+          + bcSens_[cpI].x();
         derivatives_[3*cpI + 1] =
             flowSens_[cpI].y()
           + dSdbSens_[cpI].y()
           + dndbSens_[cpI].y()
-          + dxdbDirectSens_[cpI].y();
+          + dxdbDirectSens_[cpI].y()
+          + bcSens_[cpI].y();
         derivatives_[3*cpI + 2] =
             flowSens_[cpI].z()
           + dSdbSens_[cpI].z()
           + dndbSens_[cpI].z()
-          + dxdbDirectSens_[cpI].z();
+          + dxdbDirectSens_[cpI].z()
+          + bcSens_[cpI].z();
     }
 }
 
@@ -228,6 +259,7 @@ void sensitivityVolBSplines::clearSensitivities()
     dSdbSens_ = vector::zero;
     dndbSens_ = vector::zero;
     dxdbDirectSens_ = vector::zero;
+    bcSens_ = vector::zero;
 
     SIBase::clearSensitivities();
 }
@@ -262,7 +294,10 @@ void sensitivityVolBSplines::write(const word& baseName)
             << setw(width) << "dndb::z" << " "
             << setw(width) << "dxdbDirect::x" << " "
             << setw(width) << "dxdbDirect::y" << " "
-            << setw(width) << "dxdbDirect::z" << endl;
+            << setw(width) << "dxdbDirect::z" << " "
+            << setw(width) << "dvdb::x" << " "
+            << setw(width) << "dvdb::y" << " "
+            << setw(width) << "dvdb::z" << endl;
 
         label passedCPs(0);
         label lastActive(-1);
@@ -298,7 +333,10 @@ void sensitivityVolBSplines::write(const word& baseName)
                        << setw(width) << dndbSens_[globalCP].z() << " "
                        << setw(width) << dxdbDirectSens_[globalCP].x() << " "
                        << setw(width) << dxdbDirectSens_[globalCP].y() << " "
-                       << setw(width) << dxdbDirectSens_[globalCP].z()
+                       << setw(width) << dxdbDirectSens_[globalCP].z() << " "
+                       << setw(width) << bcSens_[globalCP].x() << " "
+                       << setw(width) << bcSens_[globalCP].y() << " "
+                       << setw(width) << bcSens_[globalCP].z()
                        << endl;
                 }
             }
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplines/sensitivityVolBSplinesIncompressible.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplines/sensitivityVolBSplinesIncompressible.H
index c42dd72eb0a9ef6ed3e2e45fda5ac9dfedfe4824..5d28be24b595fa5d365287092d25968391c2b725 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplines/sensitivityVolBSplinesIncompressible.H
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplines/sensitivityVolBSplinesIncompressible.H
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -79,12 +79,16 @@ protected:
         //- on x
         vectorField dxdbDirectSens_;
 
+        //- Term dependng on the differentiation of boundary conditions
+        vectorField bcSens_;
+
         fileName derivativesFolder_;
 
 
     // Protected Member Functions
 
         void computeObjectiveContributions();
+        void computeBCContributions();
 
 
 private:
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.C
index 7191f5a90cfd571d48feefc39b01d3168a5ba5f6..3fc71d20601719bf5b8565d60051099e69844ce7 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.C
@@ -83,6 +83,7 @@ sensitivityVolBSplinesFI::sensitivityVolBSplinesFI
     dVdbSens_(0),
     distanceSens_(0),
     optionsSens_(0),
+    bcSens_(0),
 
     derivativesFolder_("optimisation"/type() + "Derivatives")
 {
@@ -97,6 +98,7 @@ sensitivityVolBSplinesFI::sensitivityVolBSplinesFI
     dVdbSens_ = vectorField(nCPs, Zero);
     distanceSens_ = vectorField(nCPs, Zero);
     optionsSens_ = vectorField(nCPs, Zero);
+    bcSens_ = vectorField(nCPs, Zero);
 
     // Create folder to store sensitivities
     mkDir(derivativesFolder_);
@@ -146,6 +148,12 @@ void sensitivityVolBSplinesFI::assembleSensitivities()
             sensitivityPatchIDs_.toc()
         );
 
+        vectorField bcSens = boxes[iNURB].computeControlPointSensitivities
+        (
+            bcDxDbMult_(),
+            sensitivityPatchIDs_.toc()
+        );
+
         for (label cpI = 0; cpI < nb; cpI++)
         {
             label globalCP = passedCPs + cpI;
@@ -267,6 +275,9 @@ void sensitivityVolBSplinesFI::assembleSensitivities()
             // dxdbSens storage
             dxdbDirectSens_[globalCP] = dxdbSens[cpI];
 
+            // bcSens storage
+            bcSens_[globalCP] = bcSens[cpI];
+
             boxSensitivities[cpI] =
                 flowSens_[globalCP]
               + dSdbSens_[globalCP]
@@ -274,7 +285,8 @@ void sensitivityVolBSplinesFI::assembleSensitivities()
               + dVdbSens_[globalCP]
               + distanceSens_[globalCP]
               + dxdbDirectSens_[globalCP]
-              + optionsSens_[globalCP];
+              + optionsSens_[globalCP]
+              + bcSens_[globalCP];
         }
 
         // Zero sensitivities in non-active design variables
@@ -302,6 +314,7 @@ void sensitivityVolBSplinesFI::assembleSensitivities()
     volBSplinesBase_.boundControlPointMovement(distanceSens_);
     volBSplinesBase_.boundControlPointMovement(dxdbDirectSens_);
     volBSplinesBase_.boundControlPointMovement(optionsSens_);
+    volBSplinesBase_.boundControlPointMovement(bcSens_);
 }
 
 
@@ -314,6 +327,7 @@ void sensitivityVolBSplinesFI::clearSensitivities()
     dVdbSens_ = vector::zero;
     distanceSens_ = vector::zero;
     optionsSens_ = vector::zero;
+    bcSens_ = vector::zero;
 
     FIBase::clearSensitivities();
 }
@@ -359,7 +373,10 @@ void sensitivityVolBSplinesFI::write(const word& baseName)
             << setw(width) << "distance::z" << " "
             << setw(width) << "options::x" << " "
             << setw(width) << "options::y" << " "
-            << setw(width) << "options::z" << endl;
+            << setw(width) << "options::z" << " "
+            << setw(width) << "dvdb::x" << " "
+            << setw(width) << "dvdb::y" << " "
+            << setw(width) << "dvdb::z" << endl;
 
         label passedCPs(0);
         label lastActive(-1);
@@ -401,7 +418,10 @@ void sensitivityVolBSplinesFI::write(const word& baseName)
                         << setw(width) << distanceSens_[globalCP].z() << " "
                         << setw(width) << optionsSens_[globalCP].x() << " "
                         << setw(width) << optionsSens_[globalCP].y() << " "
-                        << setw(width) << optionsSens_[globalCP].z() << endl;
+                        << setw(width) << optionsSens_[globalCP].z() << " "
+                        << setw(width) << bcSens_[globalCP].x() << " "
+                        << setw(width) << bcSens_[globalCP].y() << " "
+                        << setw(width) << bcSens_[globalCP].z() << endl;
                 }
             }
             passedCPs += nb;
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.H
index 24267624bafd9d2024fcbe85098ea942ba425a49..b718cdfe9348dfb0dc5ab05ab213bac8cf036cee 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.H
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivityVolBSplinesFI/sensitivityVolBSplinesFIIncompressible.H
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -88,6 +88,9 @@ protected:
         //- Term depending on fvOptions
         vectorField optionsSens_;
 
+        //- Term depending on the differentiation of boundary conditions
+        vectorField bcSens_;
+
         fileName derivativesFolder_;
 
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/shapeSensitivities/shapeSensitivitiesIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/shapeSensitivities/shapeSensitivitiesIncompressible.C
new file mode 100644
index 0000000000000000000000000000000000000000..1754b986c641afa55c9fd6ec555bb6cc5b2849f5
--- /dev/null
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/shapeSensitivities/shapeSensitivitiesIncompressible.C
@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 PCOpt/NTUA
+    Copyright (C) 2020 FOSS GP
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM 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 General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "shapeSensitivitiesIncompressible.H"
+#include "adjointBoundaryConditions.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace incompressible
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(shapeSensitivities, 0);
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void shapeSensitivities::accumulateDirectSensitivityIntegrand(const scalar dt)
+{
+    // Accumulate direct sensitivities
+    PtrList<objective>& functions(objectiveManager_.getObjectiveFunctions());
+    for (const label patchI : sensitivityPatchIDs_)
+    {
+        const scalarField magSfDt(mesh_.boundary()[patchI].magSf()*dt);
+        for (objective& func : functions)
+        {
+            const scalar wei(func.weight());
+            dSfdbMult_()[patchI] += wei*func.dSdbMultiplier(patchI)*dt;
+            dnfdbMult_()[patchI] += wei*func.dndbMultiplier(patchI)*magSfDt;
+            dxdbDirectMult_()[patchI] +=
+                wei*func.dxdbDirectMultiplier(patchI)*magSfDt;
+        }
+    }
+}
+
+
+void shapeSensitivities::accumulateBCSensitivityIntegrand(const scalar dt)
+{
+    auto& UaBoundary = adjointVars_.Ua().boundaryFieldRef();
+    tmp<boundaryVectorField> DvDbMult(dvdbMult());
+
+    // Accumulate sensitivities due to boundary conditions
+    for (const label patchI : sensitivityPatchIDs_)
+    {
+        const scalarField magSfDt(mesh_.boundary()[patchI].magSf()*dt);
+        fvPatchVectorField& Uab = UaBoundary[patchI];
+        if (isA<adjointVectorBoundaryCondition>(Uab))
+        {
+            bcDxDbMult_()[patchI] +=
+            (
+                DvDbMult()[patchI]
+              & refCast<adjointVectorBoundaryCondition>(Uab).dxdbMult()
+            )*magSfDt;
+        }
+    }
+}
+
+
+tmp<boundaryVectorField> shapeSensitivities::dvdbMult() const
+{
+    tmp<boundaryVectorField>
+        tres(createZeroBoundaryPtr<vector>(meshShape_).ptr());
+    boundaryVectorField& res = tres.ref();
+
+    // Grab references
+    const volScalarField& pa = adjointVars_.pa();
+    const volVectorField& Ua = adjointVars_.Ua();
+    const autoPtr<incompressibleAdjoint::adjointRASModel>& adjointTurbulence =
+        adjointVars_.adjointTurbulence();
+
+    // Fields needed to calculate adjoint sensitivities
+    const autoPtr<incompressible::RASModelVariables>&
+       turbVars = primalVars_.RASModelVariables();
+    const singlePhaseTransportModel& lamTransp = primalVars_.laminarTransport();
+    volScalarField nuEff(lamTransp.nu() + turbVars->nutRef());
+    volTensorField gradUa(fvc::grad(Ua));
+
+    for (const label patchI : sensitivityPatchIDs_)
+    {
+        const fvPatch& patch = meshShape_.boundary()[patchI];
+        tmp<vectorField> tnf = patch.nf();
+        const vectorField& nf = tnf();
+
+        res[patchI] =
+            (
+                nuEff.boundaryField()[patchI]
+              * (
+                    Ua.boundaryField()[patchI].snGrad()
+                  + (gradUa.boundaryField()[patchI] & nf)
+                )
+            )
+          - (nf*pa.boundaryField()[patchI])
+          + adjointTurbulence().adjointMomentumBCSource()[patchI];
+    }
+
+    return tres;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+shapeSensitivities::shapeSensitivities
+(
+    const fvMesh& mesh,
+    const dictionary& dict,
+    incompressibleVars& primalVars,
+    incompressibleAdjointVars& adjointVars,
+    objectiveManager& objectiveManager,
+    fv::optionAdjointList& fvOptionsAdjoint
+)
+:
+    adjointSensitivity
+    (
+        mesh,
+        dict,
+        primalVars,
+        adjointVars,
+        objectiveManager,
+        fvOptionsAdjoint
+    ),
+    shapeSensitivitiesBase(mesh, dict),
+    dSfdbMult_(createZeroBoundaryPtr<vector>(mesh_)),
+    dnfdbMult_(createZeroBoundaryPtr<vector>(mesh_)),
+    dxdbDirectMult_(createZeroBoundaryPtr<vector>(mesh_)),
+    bcDxDbMult_(createZeroBoundaryPtr<vector>(mesh_))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void shapeSensitivities::clearSensitivities()
+{
+    dSfdbMult_() = vector::zero;
+    dnfdbMult_() = vector::zero;
+    dxdbDirectMult_() = vector::zero;
+    bcDxDbMult_() = vector::zero;
+
+    adjointSensitivity::clearSensitivities();
+    shapeSensitivitiesBase::clearSensitivities();
+}
+
+
+void shapeSensitivities::write(const word& baseName)
+{
+    adjointSensitivity::write(baseName);
+    shapeSensitivitiesBase::write();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/shapeSensitivities/shapeSensitivitiesIncompressible.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/shapeSensitivities/shapeSensitivitiesIncompressible.H
new file mode 100644
index 0000000000000000000000000000000000000000..6b5ac239a07ac5124b6869ab222aed30be9ae451
--- /dev/null
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/shapeSensitivities/shapeSensitivitiesIncompressible.H
@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 PCOpt/NTUA
+    Copyright (C) 2020 FOSS GP
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM 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 General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::incompressible::shapeSensitivitiesBase
+
+Description
+    Base class supporting shape sensitivity derivatives for
+    incompressible flows
+
+SourceFiles
+    shapeSensitivitiesBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef shapeSensitivitiesIncompressible_H
+#define shapeSensitivitiesIncompressible_H
+
+#include "adjointSensitivityIncompressible.H"
+#include "shapeSensitivitiesBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class shapeSensitivities Declaration
+\*---------------------------------------------------------------------------*/
+
+class shapeSensitivities
+:
+    public adjointSensitivity,
+    public shapeSensitivitiesBase
+{
+protected:
+
+    // Protected data
+
+        //- Fields related to direct sensitivities
+        autoPtr<boundaryVectorField> dSfdbMult_;
+        autoPtr<boundaryVectorField> dnfdbMult_;
+        autoPtr<boundaryVectorField> dxdbDirectMult_;
+        autoPtr<boundaryVectorField> bcDxDbMult_;
+
+
+    // Protected Member Fuctions
+
+        //- Accumulate direct sensitivities
+        virtual void accumulateDirectSensitivityIntegrand(const scalar dt);
+
+        //- Accumulate sensitivities enamating from the boundary conditions
+        virtual void accumulateBCSensitivityIntegrand(const scalar dt);
+
+        //- Compute multiplier of dv_i/db
+        tmp<boundaryVectorField> dvdbMult() const;
+
+
+private:
+
+    // Private Member Functions
+
+        //- No copy construct
+        shapeSensitivities(const shapeSensitivities&) = delete;
+
+        //- No copy assignment
+        void operator=(const shapeSensitivities&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("shapeSensitivities");
+
+
+    // Constructors
+
+        //- Construct from components
+        shapeSensitivities
+        (
+            const fvMesh& mesh,
+            const dictionary& dict,
+            incompressibleVars& primalVars,
+            incompressibleAdjointVars& adjointVars,
+            objectiveManager& objectiveManager,
+            fv::optionAdjointList& fvOptionsAdjoint
+        );
+
+
+    //- Destructor
+    virtual ~shapeSensitivities() = default;
+
+
+    // Member Functions
+
+        //- Accumulate sensitivity integrands
+        virtual void accumulateIntegrand(const scalar dt) = 0;
+
+        //- Assemble sensitivities
+        virtual void assembleSensitivities() = 0;
+
+        //- Zero sensitivity fields and their constituents
+        virtual void clearSensitivities();
+
+        //- Write sensitivity fields.
+        virtual void write(const word& baseName = word::null);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace incompressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/shapeSensitivitiesBase/shapeSensitivitiesBase.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/shapeSensitivitiesBase/shapeSensitivitiesBase.C
index df928a8c319739c5f82f64f5102cdc4f9555ee56..4170963659f4f23824f08ef9e835441de4159891 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/shapeSensitivitiesBase/shapeSensitivitiesBase.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/shapeSensitivitiesBase/shapeSensitivitiesBase.C
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -172,7 +172,7 @@ void Foam::shapeSensitivitiesBase::setSensitivityPatchIDs
 }
 
 
-void Foam::shapeSensitivitiesBase::clear()
+void Foam::shapeSensitivitiesBase::clearSensitivities()
 {
     // Face-based boundary sens
     if (wallFaceSensVecPtr_.valid())
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/shapeSensitivitiesBase/shapeSensitivitiesBase.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/shapeSensitivitiesBase/shapeSensitivitiesBase.H
index 170683ccb5d5d8924c75150a75352129ef48d437..5389e6de82a74f9ac279dcd8d2cb5590400a1b60 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/shapeSensitivitiesBase/shapeSensitivitiesBase.H
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/shapeSensitivitiesBase/shapeSensitivitiesBase.H
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -179,7 +179,7 @@ public:
         void setSensitivityPatchIDs(const labelHashSet& sensPatchIDs);
 
         //- Zero sensitivity fields and their constituents
-        void clear();
+        void clearSensitivities();
 
         //- Write sensitivity fields.
         //  If valid, copies boundaryFields to volFields and writes them.