diff --git a/sdpm/src/sdpmBody.cpp b/sdpm/src/sdpmBody.cpp
index 824e89a5b0253761f0d9f70478aca39a8b4f4189..494eea5f0d4db7777e41ccc1608cd37e4affeb6f 100644
--- a/sdpm/src/sdpmBody.cpp
+++ b/sdpm/src/sdpmBody.cpp
@@ -157,6 +157,8 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
     _nLoads.resize(_nodes.size(), sdpmVector3d::Zero());
     _nLoads1Re.resize(nM, std::vector<std::vector<sdpmVector3d>>(nF, std::vector<sdpmVector3d>(_nodes.size(), sdpmVector3d::Zero())));
     _nLoads1Im.resize(nM, std::vector<std::vector<sdpmVector3d>>(nF, std::vector<sdpmVector3d>(_nodes.size(), sdpmVector3d::Zero())));
+    _iVelocityS.resize(nM, std::vector<sdpmVector3d>(getElements().size(), sdpmVector3d::Zero()));
+    _iVelocityH.resize(nM, std::vector<sdpmVector3d>(getElements().size(), sdpmVector3d::Zero()));
 }
 
 Body::~Body()
@@ -233,8 +235,8 @@ void Body::computeVelocity(size_t m, sdpmVector3d const &ui)
         ucS[i] = _motions[m]->computeSteady(*elems[i], ui);
         ucH[i] = _motions[m]->computeHarmonic(*elems[i]);
     }
-    _iVelocityS.push_back(ucS);
-    _iVelocityH.push_back(ucH);
+    _iVelocityS[m] = ucS;
+    _iVelocityH[m] = ucH;
 }
 
 /**