diff --git a/applications/solvers/stressAnalysis/solidDisplacementFoam/Make/options b/applications/solvers/stressAnalysis/solidDisplacementFoam/Make/options
index 5a13b88292..3a6868f80a 100644
--- a/applications/solvers/stressAnalysis/solidDisplacementFoam/Make/options
+++ b/applications/solvers/stressAnalysis/solidDisplacementFoam/Make/options
@@ -6,4 +6,5 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfvOptions \
     -lmeshTools
diff --git a/applications/solvers/stressAnalysis/solidDisplacementFoam/createFields.H b/applications/solvers/stressAnalysis/solidDisplacementFoam/createFields.H
index fd97d416f9..f5adee051c 100644
--- a/applications/solvers/stressAnalysis/solidDisplacementFoam/createFields.H
+++ b/applications/solvers/stressAnalysis/solidDisplacementFoam/createFields.H
@@ -1,6 +1,6 @@
 #include "readMechanicalProperties.H"
 #include "readThermalProperties.H"
-
+mesh.Sf();
 Info<< "Reading field D\n" << endl;
 volVectorField D
 (
@@ -77,3 +77,5 @@ else
 }
 
 mesh.setFluxRequired(D.name());
+
+#include "createFvOptions.H"
diff --git a/applications/solvers/stressAnalysis/solidDisplacementFoam/solidDisplacementFoam.C b/applications/solvers/stressAnalysis/solidDisplacementFoam/solidDisplacementFoam.C
index 4fa72dda9c..ae37b99386 100644
--- a/applications/solvers/stressAnalysis/solidDisplacementFoam/solidDisplacementFoam.C
+++ b/applications/solvers/stressAnalysis/solidDisplacementFoam/solidDisplacementFoam.C
@@ -36,6 +36,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "fvOptions.H"
 #include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -68,10 +69,16 @@ int main(int argc, char *argv[])
             if (thermalStress)
             {
                 volScalarField& T = Tptr();
-                solve
+                fvScalarMatrix TEqn
                 (
-                    fvm::ddt(T) == fvm::laplacian(DT, T)
+                    fvm::ddt(T) == fvm::laplacian(DT, T) + fvOptions(T)
                 );
+
+                fvOptions.constrain(TEqn);
+
+                TEqn.solve();
+
+                fvOptions.correct(T);
             }
 
             {
@@ -81,6 +88,7 @@ int main(int argc, char *argv[])
                  ==
                     fvm::laplacian(2*mu + lambda, D, "laplacian(DD,D)")
                   + divSigmaExp
+                  + fvOptions.d2dt2(D)
                 );
 
                 if (thermalStress)
@@ -89,8 +97,7 @@ int main(int argc, char *argv[])
                     DEqn += fvc::grad(threeKalpha*T);
                 }
 
-                // DEqn.setComponentReference(1, 0, vector::X, 0);
-                // DEqn.setComponentReference(1, 0, vector::Z, 0);
+                fvOptions.constrain(DEqn);
 
                 initialResidual = DEqn.solve().max().initialResidual();