diff --git a/tbox/_src/tboxw.i b/tbox/_src/tboxw.i
index 94acd7f8e464d250ca87849e88cd8cb55004a02d..c058293a7c3e98b0064ca54b5b59e80836a59063 100644
--- a/tbox/_src/tboxw.i
+++ b/tbox/_src/tboxw.i
@@ -142,8 +142,6 @@ namespace std {
 %include "wGroups.h"
 %warnfilter(+509);
 %shared_ptr(tbox::MshDeform);
-%immutable tbox::MshDeform::msh; // read-only variable (no setter)
-%immutable tbox::MshDeform::nthreads;
 %include "wMshDeform.h" // included after std templates because it needs them
 %shared_ptr(tbox::MshCrack);
 %include "wMshCrack.h"
diff --git a/tbox/src/wMshDeform.cpp b/tbox/src/wMshDeform.cpp
index bb740c91b662f7693d089fa5fae9ddb349cb1386..9e06c9c62e1ab11ded9b076c5b544c8e7aaf8e09 100644
--- a/tbox/src/wMshDeform.cpp
+++ b/tbox/src/wMshDeform.cpp
@@ -16,6 +16,7 @@
 
 #include "wMshDeform.h"
 #include "wMshData.h"
+#include "wLinearSolver.h"
 #include "wGroup.h"
 #include "wGroups.h"
 #include "wTag.h"
@@ -23,19 +24,17 @@
 #include "wElement.h"
 #include "wCache.h"
 #include "wGauss.h"
-#include "wGmres.h"
 
 #include <tbb/global_control.h>
-#include <tbb/parallel_for_each.h>
-#include <tbb/spin_mutex.h>
 #include <deque>
 
 using namespace tbox;
 
 MshDeform::MshDeform(std::shared_ptr<MshData> _msh,
-                     int _nDim, int nthrds) : wSharedObject(), nDim(_nDim),
-                                              field(false), fixed(false), moving(false),
-                                              msh(_msh), nthreads(nthrds)
+                     std::shared_ptr<tbox::LinearSolver> _linsol,
+                     int _nDim, int nthrds) : wSharedObject(), msh(_msh), linsol(_linsol),
+                                              nDim(_nDim), nthreads(nthrds),
+                                              field(false), fixed(false), moving(false)
 {
     // Check problem dimension
     if (nDim != 2 && nDim != 3)
@@ -301,14 +300,12 @@ Eigen::MatrixXd MshDeform::buildK(tbox::Element const &e, Eigen::MatrixXd const
  */
 void MshDeform::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &K)
 {
-    // Multithread
-    tbb::spin_mutex mutex;
-
     // List of triplets to build stifness matrix
     std::deque<Eigen::Triplet<double>> T;
 
     // Build stiffness matrix
-    tbb::parallel_for_each(fldElems.begin(), fldElems.end(), [&](Element *e) {
+    for (Element const *e : fldElems)
+    {
         // Hooke tensor
         Eigen::MatrixXd H;
         if (nDim == 2)
@@ -345,7 +342,6 @@ void MshDeform::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &K)
         // Elementary stiffness matrix
         Eigen::MatrixXd Ke = this->buildK(*e, H);
         // Assembly
-        tbb::spin_mutex::scoped_lock lock(mutex);
         for (size_t i = 0; i < e->nodes.size(); ++i)
         {
             for (int m = 0; m < nDim; ++m)
@@ -361,7 +357,7 @@ void MshDeform::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &K)
                 }
             }
         }
-    });
+    }
     // Periodic boundary condition step 2
     // -> for each node pair(a, b): write "pos_a - pos_b" = 0 on row a
     for (auto p : intBndPair)
@@ -446,7 +442,7 @@ void MshDeform::build(Eigen::VectorXd &f)
  */
 void MshDeform::deform()
 {
-    // Multithread
+    // Multithread (required for linear solver)
     tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
 
     // Stiffness matrix
@@ -462,18 +458,9 @@ void MshDeform::deform()
     std::cout << "done" << std::endl;
     // Solve the SoE
     std ::cout << "solving equations... " << std::flush;
-    Gmres gmres;
-    // set the preconditioner
-    gmres.setFillFactor(2);
-    gmres.setDropTol(1e-6);
-    // set the gmres
-    gmres.setTolerance(1e-8);
-    gmres.setRestart(50);
-    gmres.setMaxIterations(500);
-    // solve
     Eigen::Map<Eigen::VectorXd> f_(f.data(), f.size()), q_(q.data(), q.size());
-    gmres.compute(K, f_, q_);
-    std::cout << "done (#it: " << gmres.getIterations() << ", error: " << std::scientific << gmres.getError() << std::fixed << ")" << std::endl;
+    linsol->compute(K, f_, q_);
+    std::cout << "done (#it: " << linsol->getIterations() << ", error: " << std::scientific << linsol->getError() << std::fixed << ")" << std::endl;
 
     // Update position
     std ::cout << "updating mesh... " << std::flush;
diff --git a/tbox/src/wMshDeform.h b/tbox/src/wMshDeform.h
index e968058397e19a090196e78dadbe758523a59ad0..8f4a090c5c9a40d476be40a75f3a71ec67cd089a 100644
--- a/tbox/src/wMshDeform.h
+++ b/tbox/src/wMshDeform.h
@@ -34,8 +34,11 @@ namespace tbox
 class TBOX_API MshDeform : public fwk::wSharedObject
 {
 private:
+    std::shared_ptr<MshData> msh;                      ///< mesh
+    std::shared_ptr<tbox::LinearSolver> linsol;        ///< linear (inner) solver
     size_t mshSize;                                    ///< number of nodes in the mesh
     int nDim;                                          ///< dimension of the problem (2 or 3)
+    int nthreads;                                      ///< number of threads
     std::vector<int> rows;                             ///< unknown local index
     bool field;                                        ///< flag to check if field has been added
     bool fixed;                                        ///< flag to check if at least one fixed boundary has been added
@@ -60,10 +63,7 @@ private:
     Eigen::MatrixXd buildK(tbox::Element const &e, Eigen::MatrixXd const &H);
 
 public:
-    std::shared_ptr<MshData> msh;
-    int nthreads;
-
-    MshDeform(std::shared_ptr<MshData> _msh, int _nDim, int nthrds = 1);
+    MshDeform(std::shared_ptr<MshData> _msh, std::shared_ptr<tbox::LinearSolver> _linsol, int _nDim, int nthrds = 1);
     virtual ~MshDeform() { std::cout << "~MshDeform()\n"; }
 
     void setField(std::string const &fldName);
diff --git a/tbox/src/wUnitTest.cpp b/tbox/src/wUnitTest.cpp
index 7b886c29e77e50893e9d5c64991e1b34e538b72d..4aa67944c3db9497d80c830ba359f44918ec828d 100644
--- a/tbox/src/wUnitTest.cpp
+++ b/tbox/src/wUnitTest.cpp
@@ -78,10 +78,10 @@ void UnitTest::gmshIO(GmshExport *writer, GmshImport *reader)
 /**
  * @brief Deform a mesh of dimension d
  */
-void UnitTest::mshDeform(MshDeform *mshDef, int d)
+void UnitTest::mshDeform(MshData *msh, MshDeform *mshDef, int d)
 {
     // Initialize memory (volume element only)
-    for (auto e : mshDef->msh->elems)
+    for (auto e : msh->elems)
         if ((d == 2 && e->type() == ELTYPE::TRI3) || (d == 3 && e->type() == ELTYPE::TETRA4))
             e->initValues(true);
     // Call deform
diff --git a/tbox/src/wUnitTest.h b/tbox/src/wUnitTest.h
index 4a51567687dfde2e93d2978b576f44281b86c2ca..1170fd8e832fb573077a11034b43ebd99b48c5e0 100644
--- a/tbox/src/wUnitTest.h
+++ b/tbox/src/wUnitTest.h
@@ -36,7 +36,7 @@ public:
     virtual ~UnitTest() {}
 
     void gmshIO(GmshExport *writer, GmshImport *reader);
-    void mshDeform(MshDeform *mshDef, int d);
+    void mshDeform(MshData *msh, MshDeform *mshDef, int d);
 };
 
 } // namespace tbox
diff --git a/tbox/tests/meshDeformation.py b/tbox/tests/meshDeformation.py
index bfd2010f4832468f3b44e9bdc7b69be0a025d0bd..3de4d257b3c1256d229783ed0305725ca0c1d1ee 100644
--- a/tbox/tests/meshDeformation.py
+++ b/tbox/tests/meshDeformation.py
@@ -39,7 +39,7 @@ def main():
     msh = gmsh.MeshLoader("models/beam.geo",__file__).execute(**pars)
     gmshWriter = tbox.GmshExport(msh)
     # create mesh deformation handler
-    mshDef = tbox.MshDeform(msh, 2, nthrds=parseargs().k)
+    mshDef = tbox.MshDeform(msh, tbox.SparseLu(), 2, nthrds=parseargs().k)
     mshDef.setField("internalField")
     mshDef.addFixed(["clamp"])
     mshDef.addMoving(["tip"])
@@ -55,7 +55,7 @@ def main():
     # deform the mesh
     test = tbox.UnitTest()
     tms["deform"].start()
-    test.mshDeform(mshDef, 2)
+    test.mshDeform(msh, mshDef, 2)
     tms["deform"].stop()
     gmshWriter.save("beam_def")
 
diff --git a/tbox/tests/meshDeformation3.py b/tbox/tests/meshDeformation3.py
index 7c0599d019ff8bd58e94602b6361368f66b0face..cebff8c5ee5a6975b291be4e684d0926bd5657f3 100644
--- a/tbox/tests/meshDeformation3.py
+++ b/tbox/tests/meshDeformation3.py
@@ -38,7 +38,7 @@ def main():
     msh = gmsh.MeshLoader("models/beam3.geo",__file__).execute(**pars)
     gmshWriter = tbox.GmshExport(msh)
     # create mesh deformation handler
-    mshDef = tbox.MshDeform(msh, 3, nthrds=parseargs().k)
+    mshDef = tbox.MshDeform(msh, tbox.SparseLu(), 3, nthrds=parseargs().k)
     mshDef.setField("field")
     mshDef.addFixed(["clamp"])
     mshDef.addMoving(["tip"])
@@ -54,7 +54,7 @@ def main():
     # deform the mesh
     test = tbox.UnitTest()
     tms["deform"].start()
-    test.mshDeform(mshDef, 3)
+    test.mshDeform(msh, mshDef, 3)
     tms["deform"].stop()
     gmshWriter.save(msh.name+"_def")