diff --git a/tbox/src/wGroup.cpp b/tbox/src/wGroup.cpp
index c26634fd859ab247d3e102ea628a242f8ff84203..b271543f2441c8da0d23f96aea35a90753eaa434 100644
--- a/tbox/src/wGroup.cpp
+++ b/tbox/src/wGroup.cpp
@@ -28,7 +28,7 @@ Group::Group(std::shared_ptr<MshData> _msh, int no) : wSharedObject(), msh(_msh)
     if (it == msh->ptags.end())
     {
         std::stringstream out;
-        out << "Physical Group #" << no << " not found among physical tags";
+        out << "Group: Physical Group #" << no << " not found among physical tags";
         throw std::out_of_range(out.str());
     }
     tag = it->second;
@@ -40,8 +40,8 @@ Group::Group(std::shared_ptr<MshData> _msh, std::string const &name) : wSharedOb
     if (it == msh->ntags.end())
     {
         std::stringstream out;
-        out << "Physical Group \"" << name << "\" not found among physical tags\n";
-        out << "[info] available tags:\n\t";
+        out << "Group: Physical Group \"" << name << "\" not found among physical tags\n";
+        out << "[info] available tags: ";
         for (auto it = msh->ntags.begin(); it != msh->ntags.end(); ++it)
         {
             auto itn(it);
diff --git a/tbox/src/wMshCrack.cpp b/tbox/src/wMshCrack.cpp
index 9e11b8a8f62a14dbb8d5b511bc299795a40d08bd..60f82dfc248a140d1912b98f93471802aadbe7bf 100644
--- a/tbox/src/wMshCrack.cpp
+++ b/tbox/src/wMshCrack.cpp
@@ -28,7 +28,7 @@ using namespace tbox;
 
 MshCrack::MshCrack(std::shared_ptr<MshData> _msh, int _nDim) : wSharedObject(), msh(_msh), nDim(_nDim)
 {
-    hasCrack = true; // so that we prevent automatic data structure creation if addCrack() is not called
+    hasCrack = true; // so that we prevent automatic data structure creation if setCrack() is not called
     crckGrp = nullptr;
     exclGrp = nullptr;
 }
@@ -37,58 +37,84 @@ MshCrack::~MshCrack()
     // Delete temporary groups
     delete crckGrp;
     delete exclGrp;
-    for (size_t i = 0; i < bndGrps.size(); ++i)
-        delete bndGrps[i];
-    for (size_t i = 0; i < bndGrp.size(); ++i)
-        delete bndGrp[i];
+    for (size_t i = 0; i < sBndGrps.size(); ++i)
+        delete sBndGrps[i];
+    for (size_t i = 0; i < sBndGrp.size(); ++i)
+        delete sBndGrp[i];
+    for (size_t i = 0; i < vBndGrps.size(); ++i)
+        delete vBndGrps[i];
+    for (size_t i = 0; i < vBndGrp.size(); ++i)
+        delete vBndGrp[i];
     std::cout << "~MshCrack()\n";
 }
 
 /**
  * @brief Create temporary group for crack
  */
-void MshCrack::setCrack(std::string const &crckName)
+void MshCrack::setCrack(std::string const &name)
 {
     // Check if a crack needs to be made
-    auto it = msh->ntags.find(crckName + '_');
+    auto it = msh->ntags.find(name + '_');
     if (it != msh->ntags.end())
-        std::cout << "Groups " << crckName << " and " << crckName << "_ are already existing. MshCrack will not be run!\n";
+        std::cout << "Groups " << name << " and " << name << "_ are already existing. MshCrack will not be run\n";
     else
     {
         hasCrack = false;
-        crckGrp = new Group(msh, crckName);
+        crckGrp = new Group(msh, name);
+        if (crckGrp->tag->dim != nDim - 1)
+        {
+            std::stringstream err;
+            err << "MshCrack::setCrack: cannot create a crack of dimension " << crckGrp->tag->dim << ", dimension must be " << nDim - 1 << "!\n";
+            throw std::runtime_error(err.str());
+        }
     }
 }
 
 /**
  * @brief Create temporary group for excluded nodes on crack
  */
-void MshCrack::setExcluded(std::string const &exclName)
+void MshCrack::setExcluded(std::string const &name)
 {
     if (!hasCrack)
-        exclGrp = new Group(msh, exclName);
+        exclGrp = new Group(msh, name);
 }
 
 /**
  * @brief Create temporary groups for crack boundaries
  */
-void MshCrack::addBoundaries(std::vector<std::string> const &bndName)
+void MshCrack::addBoundary(std::string const &name)
 {
     if (!hasCrack)
     {
-        for (auto name : bndName)
+        // Check if the boundaries on both sides of the crack are explicitely provided in two different tags
+        // if so, store them
+        try
         {
-            // Check if the boundaries are explicitely separated...
-            try
+            Groups *gs = new Groups(msh, std::vector<std::string>{name, name + '_'});
+            if (gs->groups[0]->tag->dim == nDim)
+                vBndGrps.push_back(gs);
+            else if (gs->groups[0]->tag->dim == nDim - 1)
+                sBndGrps.push_back(gs);
+            else
             {
-                Groups *gs = new Groups(msh, std::vector<std::string>{name, name + '_'});
-                bndGrps.push_back(gs);
+                std::stringstream err;
+                err << "MshCrack::addBoundary: cannot add a boundary of dimension " << gs->groups[0]->tag->dim << ", dimension must be either " << nDim << " or " << nDim - 1 << "!\n";
+                throw std::runtime_error(err.str());
             }
-            // ... if not, an algorithm will later try to to separate the elements according to the side of the crack they are touching
-            catch (const std::out_of_range &)
+        }
+        // else, the boundaries are contained in a single tag, an algorithm will later try to separate the elements according to the side of the crack they are on
+        catch (const std::out_of_range &)
+        {
+            Group *g = new Group(msh, name);
+            if (g->tag->dim == nDim)
+                vBndGrp.push_back(g);
+            else if (g->tag->dim == nDim - 1)
+                sBndGrp.push_back(g);
+            else
             {
-                Group *g = new Group(msh, name);
-                bndGrp.push_back(g);
+                std::stringstream err;
+                err << "MshCrack::addBoundary: cannot add a boundary of dimension " << g->tag->dim << ", dimension must be either " << nDim << " or " << nDim - 1 << "!\n";
+                throw std::runtime_error(err.str());
             }
         }
     }
@@ -100,15 +126,14 @@ void MshCrack::addBoundaries(std::vector<std::string> const &bndName)
 void MshCrack::run()
 {
     if (hasCrack)
-        std::cout << "MshCrack not run!\n";
+        std::cout << "MshCrack not run\n";
     else
     {
         // Check that crack is bounded
-        if (bndGrps.empty() && bndGrp.empty())
-            throw std::runtime_error("MshCrack::run crack has no boundary (use MshCrack::addBoundary()!\n");
+        if ((sBndGrps.empty() && sBndGrp.empty()) || (vBndGrps.empty() && vBndGrp.empty()))
+            throw std::runtime_error("MshCrack::run, crack has no (surface and/or volume) boundary, use MshCrack::addBoundary()!\n");
         // Open the crack and modify data structure
         openCrack();
-        modifyBoundaries_();
         modifyBoundaries();
         std::cout << "Crack added on " << *(crckGrp->tag) << "(" << crckGrp->tag->name << "_ created)" << std::endl;
     }
@@ -184,7 +209,7 @@ void MshCrack::openCrack()
         else
         {
             std::stringstream err;
-            err << "Msh::Crack: could not create element of type " << e->type() << " for dimension " << nDim << "!\n";
+            err << "MshCrack::run: could not create element of type " << e->type() << " for dimension " << nDim << "!\n";
             throw std::runtime_error(err.str());
         }
         i++;
@@ -192,15 +217,123 @@ void MshCrack::openCrack()
 }
 
 /**
- * @brief Swap the elements nodes on one side of the crack and merge physical tags
+ * @brief Modify the elements on the negative side of the crack
+ *
+ * 1) Identify boundary surface elements on negative side of crack, and store their nodes
+ *   - easy case: elements are contained in different tags
+ *   - hard case: elements are contained in a single tag, they are identified by comparing the vector joining them to the crack to the normal vector of the crack (NOT ROBUST)
+ * 2) Identify boundary volume elements on negative side of crack
+ *   - easy case: elements are contained in different tags
+ *   - hard case: elements are contained in a single tag, they are identified by
+ *     - either checking if one of their nodes belong to the surface
+ *     - or, comparing the vector joining them to the crack to the normal vector of the crack
+ * 3) Swap the nodes of (surface and volume) elements located on the negative side of the crack
+ * 4) Merge the tags identifying the (surface and volume) elements on the negative side into the tags identifying the positive side
  */
-void MshCrack::modifyBoundaries_()
+void MshCrack::modifyBoundaries()
 {
-    for (auto gs : bndGrps)
+    // Store crack nodes, and compute CGs and normals if sides are contained in a single tag
+    std::unordered_set<Node *> crkNods, crkNodsAll;
+    std::vector<Eigen::Vector3d> crkCG(crckGrp->tag->elems.size()), crkN(crckGrp->tag->elems.size());
+    if (!sBndGrp.empty() || !vBndGrp.empty())
+    {
+        // nodes of crack except excluded
+        for (auto p : nodMap)
+            crkNods.insert(p.first);
+        // all nodes of crack
+        for (auto e : crckGrp->tag->elems)
+            crkNodsAll.insert(e->nodes.begin(), e->nodes.end());
+        // cg and unit normal vector
+        for (size_t i = 0; i < crckGrp->tag->elems.size(); ++i)
+        {
+            crckGrp->tag->elems[i]->computeCg();
+            crkCG[i] = crckGrp->tag->elems[i]->cg;
+            crckGrp->tag->elems[i]->computeNormal();
+            crkN[i] = crckGrp->tag->elems[i]->normal;
+        }
+    }
+
+    // Identify surface elements on negative side, and store their nodes (to find volume elements if they are contained in a single tag)
+    std::vector<Element *> _surEl;
+    std::unordered_set<Node *> _surNd;
+    // Each side is contained in a different tag
+    for (auto gs : sBndGrps)
     {
-        // Modify element nodes touching the crack
         for (auto e : gs->groups[1]->tag->elems)
-            swapNodes(e);
+        {
+            _surEl.push_back(e);
+            // only keep nodes that are not on crack (since they are on positive side)
+            for (auto n : e->nodes)
+                if (!crkNodsAll.count(n))
+                    _surNd.insert(n);
+        }
+    }
+    // Both sides are contained in a single tag
+    for (auto g : sBndGrp)
+    {
+        // find elements having at least one node on the crack, and store if on negative side
+        std::unordered_set<Element *> onCrk = findTouching(g->tag->elems, crkNods);
+        for (auto e : onCrk)
+        {
+            // find mininal distance bewteen CGs
+            e->computeCg();
+            size_t idx = findClosest(crkCG, e->cg);
+            // if dot product between crack-to-volume vector and crack normal is negative, element is on negative side
+            if ((e->cg - crkCG[idx]).dot(crkN[idx]) < 0)
+            {
+                _surEl.push_back(e);
+                // only keep nodes that are not on crack (since they are on positive side)
+                for (auto n : e->nodes)
+                    if (!crkNodsAll.count(n))
+                        _surNd.insert(n);
+            }
+        }
+    }
+
+    // Identify volume elements touching surface elements on negative side
+    std::vector<Element *> _volEl;
+    // Each side is contained in a different tag
+    for (auto gs : vBndGrps)
+        for (auto e : gs->groups[1]->tag->elems)
+            _volEl.push_back(e);
+    // Both sides are contained in a single tag
+    for (auto g : vBndGrp)
+    {
+        // find elements having at least one node on the crack
+        std::unordered_set<Element *> onCrkOnly = findTouching(g->tag->elems, crkNods);
+        // find elements having at least one node on the negative side of the surface
+        std::unordered_set<Element *> onSurf = findTouching(std::vector<Element *>(onCrkOnly.begin(), onCrkOnly.end()), _surNd);
+        // remove these elements from the first set, as they are necessarily on the negative side
+        _volEl.insert(_volEl.end(), onSurf.begin(), onSurf.end());
+        for (auto e : onSurf)
+            onCrkOnly.erase(e);
+        // check if remaining elements are on the negative side
+        for (auto e : onCrkOnly)
+        {
+            // find mininal distance bewteen CGs
+            e->computeCg();
+            size_t idx = findClosest(crkCG, e->cg);
+            // if dot product between crack-to-volume vector and crack normal is negative, element is on negative side
+            if ((e->cg - crkCG[idx]).dot(crkN[idx]) < 0)
+                _volEl.push_back(e);
+        }
+    }
+
+    // Swap nodes of surface and volume elements on negative side
+    swapNodes(_surEl);
+    swapNodes(_volEl);
+    // Merge tag identifying negative side into positive side tag
+    mergeTags(sBndGrps);
+    mergeTags(vBndGrps);
+}
+
+/**
+ * @brief Merge physical tags and modify data structure accordingly
+ */
+void MshCrack::mergeTags(std::vector<Groups *> &grps)
+{
+    for (auto gs : grps)
+    {
         // Merge tags
         Tag *refTag = gs->groups[0]->tag;
         Tag *oldTag = gs->groups[1]->tag;
@@ -217,74 +350,58 @@ void MshCrack::modifyBoundaries_()
 }
 
 /**
- * @brief Detect the elements on one side of the crack and swap their nodes
+ * @brief Find elements touching a surface (ie, that have at least a node in the provided set)
  */
-void MshCrack::modifyBoundaries()
+std::unordered_set<Element *> MshCrack::findTouching(std::vector<Element *> const &elms, std::unordered_set<Node *> const &nods)
 {
-    // compute crack CGs and normals
-    std::vector<Eigen::Vector3d> crkCG(crckGrp->tag->elems.size()), crkN(crckGrp->tag->elems.size());
-    for (size_t i = 0; i < crckGrp->tag->elems.size(); ++i)
+    std::unordered_set<Element *> found;
+    for (auto e : elms)
     {
-        crckGrp->tag->elems[i]->computeCg();
-        crkCG[i] = crckGrp->tag->elems[i]->cg;
-        crckGrp->tag->elems[i]->computeNormal();
-        crkN[i] = crckGrp->tag->elems[i]->normal;
-    }
-    // find elements on the opposite side of the crack and swap their nodes
-    for (auto g : bndGrp)
-    {
-        // find elements that have at least one node on the crack and compute their CG
-        std::map<Element *, Eigen::Vector3d> onCrk;
-        for (auto e : g->tag->elems)
-        {
-            for (auto n : e->nodes)
-            {
-                if (nodMap.find(n) != nodMap.end())
-                {
-                    //Eigen::Vector3d cg(0, 0, 0);
-                    //for (auto n_ : e->nodes)
-                    //    cg += n_->pos;
-                    //cg /= e->nodes.size();
-                    //onCrk[e] = cg;
-                    e->computeCg();
-                    onCrk[e] = e->cg;
-                    break;
-                }
-            }
-        }
-        // find closest element, check side and swap nodes
-        for (auto p : onCrk)
+        for (auto n : e->nodes)
         {
-            // find mininal distance bewteen CGs
-            size_t idx = 0;
-            std::vector<Eigen::Vector3d> deltaCG(crkCG.size());
-            for (size_t i = 0; i < crkCG.size(); ++i)
+            if (nods.count(n))
             {
-                deltaCG[i] = p.second - crkCG[i];
-                if (deltaCG[i].norm() < deltaCG[idx].norm())
-                    idx = i;
+                found.insert(e);
+                break;
             }
-            // check sign of dot product
-            if (deltaCG[idx].dot(crkN[idx]) < 0)
-                swapNodes(p.first);
         }
     }
+    return found;
 }
 
 /**
- * @brief Replace nodes of elements on opposite side of crack by opposite nodes
+ * @brief Find closest element and returns its index
  */
-void MshCrack::swapNodes(Element *e)
+size_t MshCrack::findClosest(std::vector<Eigen::Vector3d> const &cga, Eigen::Vector3d const &cgb)
 {
-    for (size_t i = 0; i < e->nodes.size(); ++i)
+    size_t idx = 0;
+    std::vector<Eigen::Vector3d> delta(cga.size());
+    for (size_t i = 0; i < cga.size(); ++i)
     {
-        try
-        {
-            e->nodes[i] = nodMap.at(e->nodes[i]);
-        }
-        catch (const std::out_of_range &)
+        delta[i] = cgb - cga[i];
+        if (delta[i].norm() < delta[idx].norm())
+            idx = i;
+    }
+    return idx;
+}
+
+/**
+ * @brief Replace nodes of boundary elements on negative side of crack by negative crack nodes
+ */
+void MshCrack::swapNodes(std::vector<Element *> &elms)
+{
+    for (auto e : elms)
+    {
+        for (size_t i = 0; i < e->nodes.size(); ++i)
         {
-            //std::cout << e->nodes[i]->no << "not found in map!\n";
+            try
+            {
+                e->nodes[i] = nodMap.at(e->nodes[i]);
+            }
+            catch (const std::out_of_range &)
+            {
+                // std::cout << e->nodes[i]->no << "not found in map!\n";
+            }
         }
     }
 }
@@ -292,4 +409,4 @@ void MshCrack::swapNodes(Element *e)
 void MshCrack::write(std::ostream &out) const
 {
     out << "MshCrack on " << msh << std::endl;
-}
\ No newline at end of file
+}
diff --git a/tbox/src/wMshCrack.h b/tbox/src/wMshCrack.h
index dcce73d0dd8de221f3ca0fa3c3c22a34ccbb5a48..92338f7c85b1f56aa76763d35794b897b0b3894e 100644
--- a/tbox/src/wMshCrack.h
+++ b/tbox/src/wMshCrack.h
@@ -19,9 +19,10 @@
 
 #include "tbox.h"
 #include "wObject.h"
-#include "wMshExport.h"
+#include <Eigen/Dense>
 #include <vector>
 #include <map>
+#include <unordered_set>
 #include <memory>
 
 namespace tbox
@@ -30,9 +31,10 @@ namespace tbox
 /**
  * @brief Duplicate nodes and elements in a physical group and make mesh consistent
  *
- * Similar to plugin "crack" from gmsh but allows more control. Currently coded
+ * Similar to plugin "crack" from gmsh but allows more control. Currently implemented
  * only for Line2 and Tri3.
  * @authors Adrien Crovato
+ * @todo Automatic identification of boundary surface elements on negative side of crack is not robust (might fail for concave boundaries)
  */
 class TBOX_API MshCrack : public fwk::wSharedObject
 {
@@ -43,21 +45,25 @@ private:
     std::map<Node *, Node *> nodMap; ///< map between nodes on both side of crack
     Group *crckGrp;                  ///< physical group of crack
     Group *exclGrp;                  ///< physical group of nodes not to be duplicated
-    std::vector<Groups *> bndGrps;   ///< list of pair of physical groups touching crack
-    std::vector<Group *> bndGrp;     ///< list of physical groups touching crack
+    std::vector<Groups *> sBndGrps;  ///< list of pair of surface physical groups touching crack
+    std::vector<Groups *> vBndGrps;  ///< list of pair of volume physical groups touching crack
+    std::vector<Group *> sBndGrp;    ///< list of surface physical groups touching crack
+    std::vector<Group *> vBndGrp;    ///< list of volume physical groups touching crack
 
     void openCrack();
-    void modifyBoundaries_();
     void modifyBoundaries();
-    void swapNodes(Element *e);
+    std::unordered_set<Element *> findTouching(std::vector<Element *> const &elms, std::unordered_set<Node *> const &nods);
+    size_t findClosest(std::vector<Eigen::Vector3d> const &cga, Eigen::Vector3d const &cgb);
+    void mergeTags(std::vector<Groups *> &grps);
+    void swapNodes(std::vector<Element *> &elms);
 
 public:
     MshCrack(std::shared_ptr<MshData> _msh, int _nDim);
     virtual ~MshCrack();
 
-    void setCrack(std::string const &crckName);
-    void addBoundaries(std::vector<std::string> const &bndName);
-    void setExcluded(std::string const &exclName);
+    void setCrack(std::string const &name);
+    void setExcluded(std::string const &name);
+    void addBoundary(std::string const &name);
     void run();
 
 #ifndef SWIG
@@ -67,4 +73,4 @@ public:
 
 } // namespace tbox
 
-#endif //WMSHCRACK_H
\ No newline at end of file
+#endif // WMSHCRACK_H
\ No newline at end of file
diff --git a/tbox/src/wMshDeform.cpp b/tbox/src/wMshDeform.cpp
index 63ca062ddca82cc32fcfc79f4cfebb36bd84b595..075a4a50e3a6d9438657c5fbaa8e8b5b7a60621e 100644
--- a/tbox/src/wMshDeform.cpp
+++ b/tbox/src/wMshDeform.cpp
@@ -34,7 +34,6 @@ MshDeform::MshDeform(std::shared_ptr<MshData> _msh,
                      std::shared_ptr<tbox::LinearSolver> _linsol,
                      int _nDim, int nthrds, int vrb) : wSharedObject(),
                                                        nDim(_nDim), nthreads(nthrds), verbose(vrb),
-                                                       field(false), fixed(false), moving(false),
                                                        msh(_msh), linsol(_linsol)
 {
     // Check problem dimension
@@ -59,37 +58,25 @@ MshDeform::MshDeform(std::shared_ptr<MshData> _msh,
  */
 void MshDeform::setField(std::string const &fldName)
 {
-    // Create temporary field group
-    fldGrp = new Group(msh, fldName);
-    field = true;
+    fldGrp = new Group(msh, fldName); // temporary field group
 }
 
 /**
  * @brief Add fixed boundaries
  */
-void MshDeform::addFixed(std::vector<std::string> const &fxdBndName)
+void MshDeform::addFixed(std::string const &fxdBndName)
 {
-    // Create temporary boundary groups
-    for (size_t i = 0; i < fxdBndName.size(); ++i)
-    {
-        Group *g = new Group(msh, fxdBndName[i]);
-        fxdBndGrp.push_back(g);
-    }
-    fixed = true;
+    Group *g = new Group(msh, fxdBndName); // temporary boundary group
+    fxdBndGrp.push_back(g);
 }
 
 /**
  * @brief Add moving boundaries
  */
-void MshDeform::addMoving(std::vector<std::string> const &movBndName)
+void MshDeform::addMoving(std::string const &movBndName)
 {
-    // Create temporary boundary groups
-    for (size_t i = 0; i < movBndName.size(); ++i)
-    {
-        Group *g = new Group(msh, movBndName[i]);
-        movBndGrp.push_back(g);
-    }
-    moving = true;
+    Group *g = new Group(msh, movBndName); // temporary boundary group
+    movBndGrp.push_back(g);
 }
 
 /**
@@ -97,23 +84,19 @@ void MshDeform::addMoving(std::vector<std::string> const &movBndName)
  */
 void MshDeform::setSymmetry(std::string const &symBndName, size_t _blck)
 {
-    // Get blocked direction
+    // Get blocked direction and create group
     blck = _blck;
     if (blck != 0 && blck != 1 && blck != 2)
-        throw std::runtime_error("tbox::MshDeform::setSymmetry: blck should be 0, 1 or 2 (for x, y and z direction resp.)!\n");
-    // Create temporary boundary groups
-    symBndGrp = new Group(msh, symBndName);
+        throw std::runtime_error("MshDeform::setSymmetry: blck should be 0, 1 or 2 (for x, y and z direction resp.)!\n");
+    symBndGrp = new Group(msh, symBndName); // temporary field group
 }
 
 /**
  * @brief Add internal boundaries
  */
-void MshDeform::addInternal(std::vector<std::string> const &intBndName)
+void MshDeform::addInternal(std::string const &intBndName, std::string const &intBndName_)
 {
-    // Create temporary internal link
-    if (intBndName.size() != 2)
-        throw std::runtime_error("MshDeform: input intBndName must be a pair of internal boundaries!\n");
-    Groups *gs = new Groups(msh, intBndName);
+    Groups *gs = new Groups(msh, std::vector<std::string>{intBndName, intBndName_}); // temporary internal link group
     intBndGrps.push_back(gs);
 }
 
@@ -123,12 +106,12 @@ void MshDeform::addInternal(std::vector<std::string> const &intBndName)
 void MshDeform::initialize()
 {
     // Sanity checks
-    if (!field)
-        throw std::runtime_error("MshDeform: field has not been found! Use setField() to add a field.\n");
-    else if (!fixed)
-        throw std::runtime_error("MshDeform: fixed boundaries have not been found! Use addFixed() to add at least one boundary.\n");
-    else if (!moving)
-        throw std::runtime_error("MshDeform: moving boundaries have not been found! Use addMoving() to add at least one boundary.\n");
+    if (fldGrp == nullptr)
+        throw std::runtime_error("MshDeform::initialize: field has not been found! Use setField() to add a field.\n");
+    if (fxdBndGrp.empty())
+        throw std::runtime_error("MshDeform::initialize: fixed boundaries have not been found! Use addFixed() to add at least one boundary.\n");
+    if (movBndGrp.empty())
+        throw std::runtime_error("MshDeform::initialize: moving boundaries have not been found! Use addMoving() to add at least one boundary.\n");
 
     // Store the field elements
     for (auto e : fldGrp->tag->elems)
@@ -186,7 +169,7 @@ void MshDeform::initialize()
         uniqueSort(nodes1);
         // Create the nodes pair
         if (nodes0.size() != nodes1.size())
-            throw std::runtime_error("MshDeform: internal boudnaries must have the same number of nodes!");
+            throw std::runtime_error("MshDeform::initialize: internal boudnaries must have the same number of nodes!");
         for (auto n0 : nodes0)
         {
             for (auto n1 : nodes1)
diff --git a/tbox/src/wMshDeform.h b/tbox/src/wMshDeform.h
index dec69f68e6e7b662b81685fb551d677113b6c79b..3d08c5709556b13364a09afbe17728232aff9dcc 100644
--- a/tbox/src/wMshDeform.h
+++ b/tbox/src/wMshDeform.h
@@ -39,9 +39,6 @@ private:
     int nthreads;                                      ///< number of threads
     int verbose;                                       ///< verbosity level
     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
-    bool moving;                                       ///< flag to check if at least one moving boundary has been added
     Group *fldGrp;                                     ///< temporary group for field
     std::vector<Group *> fxdBndGrp;                    ///< temporary group for fixed boundaries
     std::vector<Group *> movBndGrp;                    ///< temporary group for moving boundaries
@@ -69,10 +66,10 @@ public:
     virtual ~MshDeform() { std::cout << "~MshDeform()\n"; }
 
     void setField(std::string const &fldName);
-    void addMoving(std::vector<std::string> const &fxdBndName);
-    void addFixed(std::vector<std::string> const &movBndName);
+    void addMoving(std::string const &fxdBndName);
+    void addFixed(std::string const &movBndName);
     void setSymmetry(std::string const &symBndName, size_t _blck);
-    void addInternal(std::vector<std::string> const &intBndName);
+    void addInternal(std::string const &intBndName, std::string const &intBndName_);
     void initialize();
 
     void savePos();
diff --git a/tbox/tests/meshDeformation.py b/tbox/tests/meshDeformation.py
index 02c0c04dc7c5eb49449c10cf6311c43925e0d714..78f8353c8fee407ef33e4a7b55153b33366fcc60 100644
--- a/tbox/tests/meshDeformation.py
+++ b/tbox/tests/meshDeformation.py
@@ -41,8 +41,8 @@ def main():
     # create mesh deformation handler
     mshDef = tbox.MshDeform(msh, tbox.SparseLu(), 2, nthrds=parseargs().k)
     mshDef.setField("internalField")
-    mshDef.addFixed(["clamp"])
-    mshDef.addMoving(["tip"])
+    mshDef.addFixed("clamp")
+    mshDef.addMoving("tip")
     mshDef.initialize()
 
     # deform the mesh (dummy translation of "tip" nodes)
diff --git a/tbox/tests/meshDeformation3.py b/tbox/tests/meshDeformation3.py
index 61534871c7061797e911c81bb742a6d7e06d46b2..08e67b6b30f5bd0a453a2785a0a0aab61e658f28 100644
--- a/tbox/tests/meshDeformation3.py
+++ b/tbox/tests/meshDeformation3.py
@@ -40,8 +40,8 @@ def main():
     # create mesh deformation handler
     mshDef = tbox.MshDeform(msh, tbox.SparseLu(), 3, nthrds=parseargs().k)
     mshDef.setField("field")
-    mshDef.addFixed(["clamp"])
-    mshDef.addMoving(["tip"])
+    mshDef.addFixed("clamp")
+    mshDef.addMoving("tip")
     mshDef.initialize()
 
     # deform the mesh (dummy translation of "tip" nodes)