diff --git a/flow/src/wAdjoint.cpp b/flow/src/wAdjoint.cpp
index 6774d1bb73a923f294fd7e37d33670468c4bb5cb..a5d9be24ef513795d82f9b6fce75f29a7433ccca 100644
--- a/flow/src/wAdjoint.cpp
+++ b/flow/src/wAdjoint.cpp
@@ -62,7 +62,7 @@ Adjoint::Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform>
     verbose = sol->verbose;
 
     // Update element gradients
-    sol->pbl->initGradMem();
+    sol->pbl->initGradElems();
 
     // Setup adjoint variables and gradients
     lamClFlo.resize(sol->pbl->msh->nodes.size(), 0.);
diff --git a/flow/src/wBlowingResidual.cpp b/flow/src/wBlowingResidual.cpp
index 09dfcf30ec0814d621557706f060a450e658788d..150924687d7b1185ae5fc47dfa36c180a70482cb 100644
--- a/flow/src/wBlowingResidual.cpp
+++ b/flow/src/wBlowingResidual.cpp
@@ -20,7 +20,6 @@
 #include "wElement.h"
 #include "wCache.h"
 #include "wGauss.h"
-#include "wMem.h"
 
 using namespace tbox;
 using namespace flow;
@@ -35,7 +34,6 @@ Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> con
     // Get pre-computed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Blowing velocity
     double vn = f.eval(e, phi, 0);
@@ -43,6 +41,6 @@ Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> con
     // Integral of (normal) transpiration velocity
     Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
     for (size_t k = 0; k < gauss.getN(); ++k)
-        b += cache.getSf(k) * vn * gauss.getW(k) * mem.getDetJ(k);
+        b += cache.getSf(k) * vn * gauss.getW(k) * e.getDetJ(k);
     return b;
 }
diff --git a/flow/src/wBoundary.cpp b/flow/src/wBoundary.cpp
index 61e3bfb1b218df2a993fc8ab662457950bdf66ee..c01595b6fccda649f06b916a7d2f4331545cff30 100644
--- a/flow/src/wBoundary.cpp
+++ b/flow/src/wBoundary.cpp
@@ -21,10 +21,6 @@
 #include "wTag.h"
 #include "wElement.h"
 #include "wNode.h"
-#include "wFct0.h"
-#include "wCache.h"
-#include "wGauss.h"
-#include "wMem.h"
 #include <unordered_set>
 #include <iomanip>
 #include <fstream>
diff --git a/flow/src/wFreestreamResidual.cpp b/flow/src/wFreestreamResidual.cpp
index afdd43aebe18a3e3501ea4328fd8b2ac47a5fb4b..b2678b4501a4f57de19cc47146fbd58976537b30 100644
--- a/flow/src/wFreestreamResidual.cpp
+++ b/flow/src/wFreestreamResidual.cpp
@@ -20,7 +20,6 @@
 #include "wElement.h"
 #include "wCache.h"
 #include "wGauss.h"
-#include "wMem.h"
 
 using namespace tbox;
 using namespace flow;
@@ -35,15 +34,14 @@ Eigen::VectorXd FreestreamResidual::build(Element const &e, std::vector<double>
     // Get pre-computed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Normal velocity
-    double vn = fs.f.eval(e, phi, 0).dot(e.normal());
+    double vn = fs.f.eval(e, phi, 0).dot(e.getNormal());
 
     // Integral of mass flux (density times normal velocity) through freestream boundary
     Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
     for (size_t k = 0; k < gauss.getN(); ++k)
-        b += cache.getSf(k) * vn * gauss.getW(k) * mem.getDetJ(k);
+        b += cache.getSf(k) * vn * gauss.getW(k) * e.getDetJ(k);
     return b;
 }
 
@@ -57,14 +55,13 @@ Eigen::VectorXd FreestreamResidual::buildGradientAoa(tbox::Element const &e, std
     // Get pre-computed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Gradient of normal velocity
-    double dvn = fs.f.evalD(e, phi, 0).dot(e.normal());
+    double dvn = fs.f.evalD(e, phi, 0).dot(e.getNormal());
 
     // Build mass flux gradient
     Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
     for (size_t k = 0; k < gauss.getN(); ++k)
-        b += cache.getSf(k) * dvn * gauss.getW(k) * mem.getDetJ(k);
+        b += cache.getSf(k) * dvn * gauss.getW(k) * e.getDetJ(k);
     return b;
 }
diff --git a/flow/src/wKuttaElement.cpp b/flow/src/wKuttaElement.cpp
index 798fd7784d9f3e1228fefdb9ca0ec474b5b83493..c155c59f5f25a0bd908d5e792c1ec35e26c4db2d 100644
--- a/flow/src/wKuttaElement.cpp
+++ b/flow/src/wKuttaElement.cpp
@@ -16,9 +16,6 @@
 
 #include "wKuttaElement.h"
 #include "wNode.h"
-#include "wMem.h"
-#include "wCache.h"
-#include "wGauss.h"
 using namespace tbox;
 using namespace flow;
 
diff --git a/flow/src/wKuttaResidual.cpp b/flow/src/wKuttaResidual.cpp
index d7adcac995d1b079033d0e95dcab1272eaae953d..ab56bd5e95f3a9bf507f6225c62ad8d63e643441 100644
--- a/flow/src/wKuttaResidual.cpp
+++ b/flow/src/wKuttaResidual.cpp
@@ -20,7 +20,6 @@
 #include "wElement.h"
 #include "wCache.h"
 #include "wGauss.h"
-#include "wMem.h"
 
 using namespace tbox;
 using namespace flow;
@@ -35,12 +34,11 @@ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<do
     // Get shape functions, Jacobian and Gauss points
     Cache &sCache = ke.surE->getVCache();
     Gauss &sGauss = sCache.getVGauss();
-    Mem &sMem = ke.surE->getVMem();
     // Get size
     size_t nRow = ke.nRow;
 
     // Velocity
-    Eigen::RowVectorXd Aup = ke.volE->computeGradient(phi, 0).transpose() * ke.volE->getVMem().getJinv(0) * ke.volE->getVCache().getDsf(0);
+    Eigen::RowVectorXd Aup = ke.volE->computeGradient(phi, 0).transpose() * ke.volE->getJinv(0) * ke.volE->getVCache().getDsf(0);
 
     // Build
     Eigen::MatrixXd A = Eigen::MatrixXd::Zero(nRow, ke.nCol);
@@ -52,7 +50,7 @@ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<do
         for (size_t i = 0; i < nRow; ++i)
             sfp(i) = sf(ke.teN[i].first);
         // wake contribution
-        A += ke.sign * sfp * Aup * sGauss.getW(k) * sMem.getDetJ(k);
+        A += ke.sign * sfp * Aup * sGauss.getW(k) * ke.surE->getDetJ(k);
     }
     return A;
 }
@@ -67,7 +65,6 @@ Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double>
     // Get shape functions, Jacobian and Gauss points
     Cache &sCache = ke.surE->getVCache();
     Gauss &sGauss = sCache.getVGauss();
-    Mem &sMem = ke.surE->getVMem();
     // Get size
     size_t nRow = ke.nRow;
 
@@ -85,7 +82,7 @@ Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double>
         for (size_t i = 0; i < nRow; ++i)
             sfp(i) = sf(ke.teN[i].first);
         // wake contribution
-        b += ke.sign * sfp * dPhi2 * sGauss.getW(k) * sMem.getDetJ(k);
+        b += ke.sign * sfp * dPhi2 * sGauss.getW(k) * ke.surE->getDetJ(k);
     }
     return b;
 }
@@ -115,10 +112,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku
     Cache &surCache = ke.surE->getVCache();
     Gauss &surGauss = surCache.getVGauss();
     // Get Jacobian
-    Mem &surMem = ke.surE->getVMem();
-    Mem &volMem = ke.volE->getVMem();
-    Eigen::MatrixXd const &iJ = volMem.getJinv(0);
-    std::vector<Eigen::MatrixXd> const &dJ = volMem.getGradJ(0);
+    Eigen::MatrixXd const &iJ = ke.volE->getJinv(0);
+    std::vector<Eigen::MatrixXd> const &dJ = ke.volE->getGradJ(0);
     // Get size
     size_t nRow = ke.nRow;
     size_t nDim = ke.nDim;
@@ -135,8 +130,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku
     {
         // weight and Jacobian determinant
         double w = surGauss.getW(k);
-        double detJ = surMem.getDetJ(k);
-        std::vector<double> dDetJ = surMem.getGradDetJ(k);
+        double detJ = ke.surE->getDetJ(k);
+        std::vector<double> dDetJ = ke.surE->getGradDetJ(k);
         // shape functions
         Eigen::VectorXd const &sf = surCache.getSf(k);
         Eigen::VectorXd sfp(nRow);
diff --git a/flow/src/wLoadFunctional.cpp b/flow/src/wLoadFunctional.cpp
index feefdaf9b2fdfc89b04b9fcb2dec2dd8a58acd17..6c0bbe17b7325badaadfe0f851fc64b87873e393 100644
--- a/flow/src/wLoadFunctional.cpp
+++ b/flow/src/wLoadFunctional.cpp
@@ -20,7 +20,6 @@
 #include "wElement.h"
 #include "wCache.h"
 #include "wGauss.h"
-#include "wMem.h"
 
 using namespace tbox;
 using namespace flow;
@@ -35,17 +34,16 @@ Eigen::VectorXd LoadFunctional::build(Element const &e, Element const &eV, std::
     // Get pre-computed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Pressure force coefficient in associated volume
-    Eigen::Vector3d fcp = cp.eval(eV, phi, 0) * e.normal();
+    Eigen::Vector3d fcp = cp.eval(eV, phi, 0) * e.getNormal();
 
     // Integrate cp on the surface
     Eigen::VectorXd b = Eigen::VectorXd::Zero(3 * e.nodes.size());
     for (size_t k = 0; k < gauss.getN(); ++k)
     {
         Eigen::VectorXd sf = cache.getSf(k);
-        double wdj = gauss.getW(k) * mem.getDetJ(k);
+        double wdj = gauss.getW(k) * e.getDetJ(k);
         for (size_t i = 0; i < e.nodes.size(); ++i)
             b.block<3, 1>(i * 3, 0) += sf(i) * fcp * wdj; // stack nodal force vectors [fx1, fy1, fz1, fx2, fy2, fz2, ...]
     }
@@ -62,17 +60,16 @@ Eigen::MatrixXd LoadFunctional::buildGradientFlow(Element const &e, Element cons
     // Get pre-computed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Gradient of pressure force coefficient in associated volume
-    Eigen::MatrixXd dfcp = e.normal() * cp.evalD(eV, phi, 0) * eV.computeGradient(phi, 0).transpose() * eV.getVMem().getJinv(0) * eV.getVCache().getDsf(0);
+    Eigen::MatrixXd dfcp = e.getNormal() * cp.evalD(eV, phi, 0) * eV.computeGradient(phi, 0).transpose() * eV.getJinv(0) * eV.getVCache().getDsf(0);
 
     // Integrate cp on the surface
     Eigen::MatrixXd A = Eigen::MatrixXd::Zero(3 * e.nodes.size(), eV.nodes.size());
     for (size_t k = 0; k < gauss.getN(); ++k)
     {
         Eigen::VectorXd sf = cache.getSf(k);
-        double wdj = gauss.getW(k) * mem.getDetJ(k);
+        double wdj = gauss.getW(k) * e.getDetJ(k);
         for (size_t i = 0; i < e.nodes.size(); ++i)
             A.block(i * 3, 0, 3, eV.nodes.size()) += sf(i) * dfcp * wdj;
     }
@@ -92,17 +89,16 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> LoadFunctional::buildGradientMesh(E
     // Get pre-computed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Variables
-    double cP = cp.eval(eV, phi, 0);  // pressure coefficient
-    Eigen::Vector3d nrm = e.normal(); // unit normal
+    double cP = cp.eval(eV, phi, 0);     // pressure coefficient
+    Eigen::Vector3d nrm = e.getNormal(); // unit normal
     // Gradients
-    double dCp = cp.evalD(eV, phi, 0);                                 // pressure coefficient
-    Eigen::VectorXd dPhi = eV.computeGradient(phi, 0);                     // potential
-    Eigen::MatrixXd const &iJ = eV.getVMem().getJinv(0);               // inverse Jacobian
-    std::vector<Eigen::MatrixXd> const &dJ = eV.getVMem().getGradJ(0); // Jacobian
-    std::vector<Eigen::Vector3d> dNrm = e.gradientNormal();            // unit normal
+    double dCp = cp.evalD(eV, phi, 0);                       // pressure coefficient
+    Eigen::VectorXd dPhi = eV.computeGradient(phi, 0);       // potential
+    Eigen::MatrixXd const &iJ = eV.getJinv(0);               // inverse Jacobian
+    std::vector<Eigen::MatrixXd> const &dJ = eV.getGradJ(0); // Jacobian
+    std::vector<Eigen::Vector3d> dNrm = e.getGradNormal();   // unit normal
 
     // Build element contributions
     Eigen::MatrixXd As = Eigen::MatrixXd::Zero(3 * e.nodes.size(), e.nodes.size() * nDim);
@@ -111,8 +107,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> LoadFunctional::buildGradientMesh(E
     {
         Eigen::VectorXd sf = cache.getSf(k);
         double w = gauss.getW(k);
-        double detJ = mem.getDetJ(k);
-        std::vector<double> dDetJ = mem.getGradDetJ(k);
+        double detJ = e.getDetJ(k);
+        std::vector<double> dDetJ = e.getGradDetJ(k);
         for (size_t i = 0; i < e.nodes.size(); ++i)
         {
             for (size_t m = 0; m < nDim; ++m)
diff --git a/flow/src/wNewton.cpp b/flow/src/wNewton.cpp
index 71a50f75ca6e842d13d9e4d35807aa3d217c4520..ec0a579a6825e21c5ae398317124b26d45443904 100644
--- a/flow/src/wNewton.cpp
+++ b/flow/src/wNewton.cpp
@@ -33,7 +33,6 @@
 #include "wTag.h"
 #include "wNode.h"
 #include "wElement.h"
-#include "wMem.h"
 #include "wLinesearch.h"
 #include "wLinearSolver.h"
 
@@ -508,10 +507,10 @@ void Newton::findUpwd()
         // -> find (cgU-cgE) \dot vel_ closest to one
         double SPbest = -1;
         size_t idx = 0;
-        Eigen::VectorXd cgE = (p.first)->getVMem().getCg().block(0, 0, pbl->nDim, 1);
+        Eigen::VectorXd cgE = (p.first)->getCg().block(0, 0, pbl->nDim, 1);
         for (size_t i = 0; i < p.second.size(); ++i)
         {
-            Eigen::VectorXd cgU = (p.second[i])->getVMem().getCg().block(0, 0, pbl->nDim, 1);
+            Eigen::VectorXd cgU = (p.second[i])->getCg().block(0, 0, pbl->nDim, 1);
             double SP = (cgU - cgE).normalized().dot(vel_);
             if (fabs(SP - 1) < fabs(SPbest - 1))
             {
diff --git a/flow/src/wPotentialResidual.cpp b/flow/src/wPotentialResidual.cpp
index f98adbc55f97dd7e3126fb12c4e7c300239e2570..1e683747272a8b1da7f599e23c1e482dd2d7f0e4 100644
--- a/flow/src/wPotentialResidual.cpp
+++ b/flow/src/wPotentialResidual.cpp
@@ -20,7 +20,6 @@
 #include "wElement.h"
 #include "wCache.h"
 #include "wGauss.h"
-#include "wMem.h"
 
 using namespace tbox;
 using namespace flow;
@@ -35,16 +34,15 @@ Eigen::MatrixXd PotentialResidual::buildFixed(Element const &e, std::vector<doub
     // Get pre-computed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Build
     Eigen::MatrixXd A = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
     for (size_t k = 0; k < gauss.getN(); ++k)
     {
         // Shape functions gradient
-        Eigen::MatrixXd const &dSf = mem.getJinv(k) * cache.getDsf(k);
+        Eigen::MatrixXd const &dSf = e.getJinv(k) * cache.getDsf(k);
         // Elementary stiffness matrix
-        A += dSf.transpose() * dSf * fluid.rho.eval(e, phi, k) * gauss.getW(k) * mem.getDetJ(k);
+        A += dSf.transpose() * dSf * fluid.rho.eval(e, phi, k) * gauss.getW(k) * e.getDetJ(k);
     }
     return A;
 }
@@ -60,12 +58,11 @@ Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, st
     // Get pre-computed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Subsonic contribution
     Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
     for (size_t k = 0; k < gauss.getN(); ++k)
-        b += fluid.rho.eval(e, phi, k) * (mem.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * mem.getDetJ(k);
+        b += fluid.rho.eval(e, phi, k) * (e.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * e.getDetJ(k);
 
     // Supersonic contribution
     double mach = fluid.mach.eval(e, phi, 0);
@@ -77,7 +74,7 @@ Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, st
         b *= 1 - mu;
         // contribution of upwind element
         for (size_t k = 0; k < gauss.getN(); ++k)
-            b += mu * rhoU * (mem.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * mem.getDetJ(k);
+            b += mu * rhoU * (e.getJinv(k) * cache.getDsf(k)).transpose() * e.computeGradient(phi, k) * gauss.getW(k) * e.getDetJ(k);
     }
     return b;
 }
@@ -97,16 +94,15 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
     // Get pre-computed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Subsonic contribution
     Eigen::MatrixXd A1 = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
     for (size_t k = 0; k < gauss.getN(); ++k)
     {
         // Gauss point and determinant
-        double wdj = gauss.getW(k) * mem.getDetJ(k);
+        double wdj = gauss.getW(k) * e.getDetJ(k);
         // Shape functions and solution gradient
-        Eigen::MatrixXd const &dSf = mem.getJinv(k) * cache.getDsf(k);
+        Eigen::MatrixXd const &dSf = e.getJinv(k) * cache.getDsf(k);
         Eigen::VectorXd dPhi = e.computeGradient(phi, k);
 
         // rho * grad_phi*grad_psi + drho * grad_phi*grad_psi
@@ -126,7 +122,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
         Eigen::VectorXd dPhiU = eU.computeGradient(phi, 0);
         double rhoU = fluid.rho.eval(eU, phi, 0);
         double dRhoU = fluid.rho.evalD(eU, phi, 0);
-        Eigen::MatrixXd const &dSfU = eU.getVMem().getJinv(0) * eU.getVCache().getDsf(0);
+        Eigen::MatrixXd const &dSfU = eU.getJinv(0) * eU.getVCache().getDsf(0);
 
         // upwinding
         A1 *= 1 - mu;
@@ -134,9 +130,9 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
         for (size_t k = 0; k < gauss.getN(); ++k)
         {
             // Gauss point and determinant
-            double wdj = gauss.getW(k) * mem.getDetJ(k);
+            double wdj = gauss.getW(k) * e.getDetJ(k);
             // Shape functions and solution gradient
-            Eigen::MatrixXd const &dSf = mem.getJinv(k) * cache.getDsf(k);
+            Eigen::MatrixXd const &dSf = e.getJinv(k) * cache.getDsf(k);
             Eigen::VectorXd dPhi = e.computeGradient(phi, k);
 
             // mu*drhoU*grad_phi*grad_psi
@@ -166,7 +162,6 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
     // Get pre-computed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Subsonic contribution
     Eigen::MatrixXd A1 = Eigen::MatrixXd::Zero(e.nodes.size(), nDim * e.nodes.size());
@@ -177,13 +172,13 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
         double rho = fluid.rho.eval(e, phi, k);
         double dRho = fluid.rho.evalD(e, phi, k);
         Eigen::MatrixXd const &dSf = cache.getDsf(k);
-        Eigen::MatrixXd const &iJ = mem.getJinv(k);
+        Eigen::MatrixXd const &iJ = e.getJinv(k);
         // Jacobian gradients
-        std::vector<Eigen::MatrixXd> const &dJ = mem.getGradJ(k);
-        std::vector<double> const &dDetJ = mem.getGradDetJ(k);
+        std::vector<Eigen::MatrixXd> const &dJ = e.getGradJ(k);
+        std::vector<double> const &dDetJ = e.getGradDetJ(k);
         // Gauss points and determinant
         double w = gauss.getW(k);
-        double detJ = mem.getDetJ(k);
+        double detJ = e.getDetJ(k);
         for (size_t i = 0; i < e.nodes.size(); ++i)
         {
             for (size_t m = 0; m < nDim; ++m)
@@ -208,8 +203,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
         Eigen::VectorXd dPhiU = eU.computeGradient(phi, 0);
         double rhoU = fluid.rho.eval(eU, phi, 0);
         double dRhoU = fluid.rho.evalD(eU, phi, 0);
-        Eigen::MatrixXd const &iJU = eU.getVMem().getJinv(0);
-        std::vector<Eigen::MatrixXd> const &dJU = eU.getVMem().getGradJ(0);
+        Eigen::MatrixXd const &iJU = eU.getJinv(0);
+        std::vector<Eigen::MatrixXd> const &dJU = eU.getGradJ(0);
 
         // upwinding
         A1 *= 1 - mu;
@@ -221,13 +216,13 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMes
             double rho = fluid.rho.eval(e, phi, k);
             double dM = fluid.mach.evalD(e, phi, k);
             Eigen::MatrixXd const &dSf = cache.getDsf(k);
-            Eigen::MatrixXd const &iJ = mem.getJinv(k);
+            Eigen::MatrixXd const &iJ = e.getJinv(k);
             // Jacobian gradients
-            std::vector<Eigen::MatrixXd> const &dJ = mem.getGradJ(k);
-            std::vector<double> const &dDetJ = mem.getGradDetJ(k);
+            std::vector<Eigen::MatrixXd> const &dJ = e.getGradJ(k);
+            std::vector<double> const &dDetJ = e.getGradDetJ(k);
             // Gauss points and determinant
             double w = gauss.getW(k);
-            double detJ = mem.getDetJ(k);
+            double detJ = e.getDetJ(k);
             for (size_t m = 0; m < nDim; ++m)
             {
                 for (size_t i = 0; i < e.nodes.size(); ++i)
diff --git a/flow/src/wProblem.cpp b/flow/src/wProblem.cpp
index 966064cab05ee1e7b50f696c12745d3d5448bd0d..7d3d069fe05aa57a3b2cc416cff3b621e52026de 100644
--- a/flow/src/wProblem.cpp
+++ b/flow/src/wProblem.cpp
@@ -16,7 +16,6 @@
 
 #include "wProblem.h"
 #include "wElement.h"
-#include "wMem.h"
 #include "wMedium.h"
 #include "wBoundary.h"
 #include "wAssign.h"
@@ -181,63 +180,62 @@ void Problem::add(std::shared_ptr<Blowing> b)
 }
 
 /**
- * @brief Initialize the elements memory (values)
- * @authors Adrien Crovato
+ * @brief Initialize the elements precomputed values
  */
-void Problem::initMem()
+void Problem::initElems()
 {
     // Update volume
     for (auto e : medium->tag->elems)
-        e->getVMem().update(true);
+        e->initValues(true);
     // Update surface (boundaries)
     for (auto surf : bnds)
         for (auto e : surf->groups[0]->tag->elems)
-            e->getVMem().update(false);
+            e->initValues(false);
     // Update surface (Freestream B.C.)
     for (auto surf : fBCs)
         for (auto e : surf->tag->elems)
-            e->getVMem().update(false);
+            e->initValues(false);
     // Update surface (Wake B.C.)
     for (auto surf : wBCs)
     {
         for (auto e : surf->groups[0]->tag->elems)
-            e->getVMem().update(false);
+            e->initValues(false);
         for (auto e : surf->groups[1]->tag->elems)
-            e->getVMem().update(false);
+            e->initValues(false);
     }
     // Update surface (Blowing B.C.) if not already done
     for (auto surf : bBCs)
         for (auto e : surf->tag->elems)
-            if (!(e->getVMem().isValid()))
-                e->getVMem().update(false);
+            if (!(e->hasValues()))
+                e->initValues(false);
 }
 
 /**
- * @brief Initialize the elements memory (gradients)
+ * @brief Initialize the elements precomputed gradients
  * @authors Adrien Crovato
  */
-void Problem::initGradMem()
+void Problem::initGradElems()
 {
     // Update volume
     for (auto e : medium->tag->elems)
-        e->getVMem().initGradients(true);
+        e->initGradients();
     // Update surface (boundaries)
     for (auto surf : bnds)
         for (auto e : surf->groups[0]->tag->elems)
-            e->getVMem().initGradients(false);
+            e->initGradients();
     // Update surface (Wake B.C.)
     for (auto surf : wBCs)
     {
         for (auto e : surf->groups[0]->tag->elems)
-            e->getVMem().initGradients(false);
+            e->initGradients();
         for (auto e : surf->groups[1]->tag->elems)
-            e->getVMem().initGradients(false);
+            e->initGradients();
     }
     // Update surface (Blowing B.C.) if not already done
     for (auto surf : bBCs)
         for (auto e : surf->tag->elems)
-            if (!(e->getVMem().isGradientValid()))
-                e->getVMem().initGradients(false);
+            if (!(e->hasGradients()))
+                e->initGradients();
 }
 
 /**
diff --git a/flow/src/wProblem.h b/flow/src/wProblem.h
index 3f4b9b630d5f130dda2744ebc9b0f8079b54596d..152aa1f267cc0530190b5c386f8ecf2b591ebac6 100644
--- a/flow/src/wProblem.h
+++ b/flow/src/wProblem.h
@@ -78,8 +78,8 @@ public:
 #ifndef SWIG
     void check() const;
     void pin();
-    void initMem();
-    void initGradMem();
+    void initElems();
+    void initGradElems();
     virtual void write(std::ostream &out) const override;
 #endif
 };
diff --git a/flow/src/wSolver.cpp b/flow/src/wSolver.cpp
index d44e29dcfdaa7af01618a2e658dfe728fec40961..1b15b0495f31c1a4316d4127ee95800d050293d5 100644
--- a/flow/src/wSolver.cpp
+++ b/flow/src/wSolver.cpp
@@ -73,7 +73,7 @@ Solver::Solver(std::shared_ptr<Problem> _pbl, std::shared_ptr<LinearSolver> _lin
     // Check the problem, pin a degree of freedom if necessary, and update element memory
     pbl->check();
     pbl->pin();
-    pbl->initMem();
+    pbl->initElems();
 
     // Assign a row (unknown index) for each node, and setup the local list
     rows.resize(pbl->msh->nodes.size());
diff --git a/flow/src/wWakeElement.cpp b/flow/src/wWakeElement.cpp
index a75ccaa0e91a73580e0dc1c3d8e2dcffc159bf38..6146ed78cd565f4df6d6e7e1f0cde6f997b0e701 100644
--- a/flow/src/wWakeElement.cpp
+++ b/flow/src/wWakeElement.cpp
@@ -16,9 +16,6 @@
 
 #include "wWakeElement.h"
 #include "wNode.h"
-#include "wMem.h"
-#include "wCache.h"
-#include "wGauss.h"
 using namespace tbox;
 using namespace flow;
 
diff --git a/flow/src/wWakeResidual.cpp b/flow/src/wWakeResidual.cpp
index 0b111e233d2e82a703972e027aa134bbe739a014..12562c7ca3b69befdb65c2c0b6e4475e8ededd96 100644
--- a/flow/src/wWakeResidual.cpp
+++ b/flow/src/wWakeResidual.cpp
@@ -20,7 +20,6 @@
 #include "wElement.h"
 #include "wCache.h"
 #include "wGauss.h"
-#include "wMem.h"
 
 using namespace tbox;
 using namespace flow;
@@ -40,9 +39,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElemen
     Eigen::MatrixXd const &vUpDsf = we.volUpE->getVCache().getDsf(0);
     Eigen::MatrixXd const &vLwDsf = we.volLwE->getVCache().getDsf(0);
     // Get Jacobian
-    Mem &sMem = we.surUpE->getVMem();
-    Eigen::MatrixXd const &vUpJ = we.volUpE->getVMem().getJinv(0);
-    Eigen::MatrixXd const &vLwJ = we.volLwE->getVMem().getJinv(0);
+    Eigen::MatrixXd const &vUpJ = we.volUpE->getJinv(0);
+    Eigen::MatrixXd const &vLwJ = we.volLwE->getJinv(0);
     // Get size
     size_t nRow = we.nRow;
     size_t nDim = we.nDim;
@@ -53,7 +51,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElemen
         for (size_t j = 0; j < nDim; ++j)
             dSf(i, j) = vUpDsf(j, we.vsMap.at(we.wkN[i].first));
     Eigen::VectorXd dSfp(nRow);
-    dSfp = 0.5 * sqrt(sMem.getVol()) * dSf * vUpJ.transpose() * vi;
+    dSfp = 0.5 * sqrt(we.surUpE->getVol()) * dSf * vUpJ.transpose() * vi;
     // Velocity
     Eigen::RowVectorXd Aup = we.volUpE->computeGradient(phi, 0).transpose() * vUpJ * vUpDsf;
     Eigen::RowVectorXd Alw = we.volLwE->computeGradient(phi, 0).transpose() * vLwJ * vLwDsf;
@@ -69,8 +67,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElemen
         for (size_t i = 0; i < nRow; ++i)
             sfp(i) = dSfp(i) + sf(we.wkN[i].first);
         // wake contribution
-        Au += sfp * Aup * sGauss.getW(k) * sMem.getDetJ(k);
-        Al += sfp * Alw * sGauss.getW(k) * sMem.getDetJ(k);
+        Au += sfp * Aup * sGauss.getW(k) * we.surUpE->getDetJ(k);
+        Al += sfp * Alw * sGauss.getW(k) * we.surUpE->getDetJ(k);
     }
     return std::make_tuple(Au, Al);
 }
@@ -87,8 +85,6 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con
     Cache &sCache = we.surUpE->getVCache();
     Gauss &sGauss = sCache.getVGauss();
     Eigen::MatrixXd const &vUpDsf = we.volUpE->getVCache().getDsf(0);
-    // Get Jacobian
-    Mem &sMem = we.surUpE->getVMem();
     // Get size
     size_t nRow = we.nRow;
     size_t nDim = we.nDim;
@@ -99,7 +95,7 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con
         for (size_t j = 0; j < nDim; ++j)
             dSf(i, j) = vUpDsf(j, we.vsMap.at(we.wkN[i].first));
     Eigen::VectorXd dSfp(nRow);
-    dSfp = 0.5 * sqrt(sMem.getVol()) * dSf * we.volUpE->getVMem().getJinv(0).transpose() * vi;
+    dSfp = 0.5 * sqrt(we.surUpE->getVol()) * dSf * we.volUpE->getJinv(0).transpose() * vi;
     // Velocity squared
     Eigen::VectorXd dPhiU = we.volUpE->computeGradient(phi, 0);
     Eigen::VectorXd dPhiL = we.volLwE->computeGradient(phi, 0);
@@ -117,8 +113,8 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con
         for (size_t i = 0; i < nRow; ++i)
             sfp(i) = dSfp(i) + sf(we.wkN[i].first);
         // wake contribution
-        bu += sfp * dPhi2u * sGauss.getW(k) * sMem.getDetJ(k);
-        bl += sfp * dPhi2l * sGauss.getW(k) * sMem.getDetJ(k);
+        bu += sfp * dPhi2u * sGauss.getW(k) * we.surUpE->getDetJ(k);
+        bl += sfp * dPhi2l * sGauss.getW(k) * we.surUpE->getDetJ(k);
     }
     return std::make_tuple(bu, bl);
 }
@@ -152,15 +148,12 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buil
     Gauss &surGauss = sCache.getVGauss();
     Eigen::MatrixXd const &dSfUp = we.volUpE->getVCache().getDsf(0);
     // Get Jacobian
-    Mem &sMem = we.surUpE->getVMem();
-    Mem &vUpMem = we.volUpE->getVMem();
-    Mem &vLwMem = we.volLwE->getVMem();
-    Eigen::MatrixXd const &iJUp = vUpMem.getJinv(0);
-    Eigen::MatrixXd const &iJLw = vLwMem.getJinv(0);
-    std::vector<Eigen::MatrixXd> const &dJUp = vUpMem.getGradJ(0);
-    std::vector<Eigen::MatrixXd> const &dJLw = vLwMem.getGradJ(0);
+    Eigen::MatrixXd const &iJUp = we.volUpE->getJinv(0);
+    Eigen::MatrixXd const &iJLw = we.volLwE->getJinv(0);
+    std::vector<Eigen::MatrixXd> const &dJUp = we.volUpE->getGradJ(0);
+    std::vector<Eigen::MatrixXd> const &dJLw = we.volLwE->getGradJ(0);
     // Get surface gradient
-    std::vector<double> const &dSurf = sMem.getGradVol();
+    std::vector<double> const &dSurf = we.surUpE->getGradVol();
     // Get size
     size_t nRow = we.nRow;
     size_t nDim = we.nDim;
@@ -169,7 +162,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buil
     size_t nSur = we.surUpE->nodes.size();
 
     // Stabilization term
-    double sqSurf = sqrt(sMem.getVol());
+    double sqSurf = sqrt(we.surUpE->getVol());
     Eigen::MatrixXd dSf(nRow, nDim);
     for (size_t i = 0; i < nRow; ++i)
         for (size_t m = 0; m < nDim; ++m)
@@ -206,8 +199,8 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buil
     {
         // weight and Jacobian determinant
         double w = surGauss.getW(k);
-        double detJ = sMem.getDetJ(k);
-        std::vector<double> dDetJ = sMem.getGradDetJ(k);
+        double detJ = we.surUpE->getDetJ(k);
+        std::vector<double> dDetJ = we.surUpE->getGradDetJ(k);
         // shape functions and stabilization
         Eigen::VectorXd const &sf = sCache.getSf(k);
         Eigen::VectorXd sfp(nRow);
diff --git a/heat/src/wHeatTerm.cpp b/heat/src/wHeatTerm.cpp
index 057e139c461a5f2b5d2b2a240d158a8ef448db9b..6d0238a6214e49f2c5aa8dbe5aae09321b2880db 100644
--- a/heat/src/wHeatTerm.cpp
+++ b/heat/src/wHeatTerm.cpp
@@ -19,7 +19,6 @@
 #include "wElement.h"
 #include "wCache.h"
 #include "wGauss.h"
-#include "wMem.h"
 #include "wFct0.h"
 #include "wFct2.h"
 
@@ -47,17 +46,16 @@ Eigen::MatrixXd HeatTerm::build(Element const &e, std::vector<double> const &u,
     else // true run - use MPI cached results
     {
         // Gauss integration
-        Mem &mem = e.getVMem();
         for (size_t k = 0; k < gauss.getN(); ++k)
         {
             // Jacobian inverse, shape functions and function evaluation
-            Eigen::MatrixXd const &J = mem.getJinv(k);
+            Eigen::MatrixXd const &J = e.getJinv(k);
             Eigen::MatrixXd const &dff = cache.getDsf(k);
             Eigen::MatrixXd fk;
             f.eval(e, u, k, fk, fake);
 
             // Elementary stiffness matrix
-            K += (fk * J * dff).transpose() * J * dff * gauss.getW(k) * mem.getDetJ(k);
+            K += (fk * J * dff).transpose() * J * dff * gauss.getW(k) * e.getDetJ(k);
         }
     }
     return K;
@@ -71,12 +69,11 @@ Eigen::VectorXd HeatTerm::build(Element const &e, std::vector<double> const &u,
     // Get precomputed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Gauss integration
     Eigen::VectorXd s = Eigen::VectorXd::Zero(e.nodes.size());
     for (size_t k = 0; k < gauss.getN(); ++k)
-        s += cache.getSf(k) * f.eval(e, u, k) * gauss.getW(k) * mem.getDetJ(k);
+        s += cache.getSf(k) * f.eval(e, u, k) * gauss.getW(k) * e.getDetJ(k);
     return s;
 }
 
@@ -89,7 +86,6 @@ Eigen::VectorXd HeatTerm::build2(Element const &e, std::vector<double> const &u,
     // Get precomputed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Gauss integration
     Eigen::VectorXd qV = Eigen::VectorXd::Zero(e.nodes.size());
@@ -101,7 +97,7 @@ Eigen::VectorXd HeatTerm::build2(Element const &e, std::vector<double> const &u,
         f.eval(e, u, k, fk, false);
 
         // Elementary flux vector
-        qV += (fk * e.computeGradient(u, k)).transpose() * mem.getJinv(k) * dff * gauss.getW(k) * mem.getDetJ(k);
+        qV += (fk * e.computeGradient(u, k)).transpose() * e.getJinv(k) * dff * gauss.getW(k) * e.getDetJ(k);
     }
     return qV;
 }
@@ -114,7 +110,6 @@ Eigen::VectorXd HeatTerm::computeFlux(Element const &e, std::vector<double> cons
     // Get precomputed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Gauss integration
     Eigen::VectorXd qV = Eigen::Vector2d::Zero(); // 2 dimensions
@@ -125,7 +120,7 @@ Eigen::VectorXd HeatTerm::computeFlux(Element const &e, std::vector<double> cons
         Eigen::MatrixXd fk(2, 2);
         f.eval(e, u, k, fk, false);
 
-        qV -= fk * gradk * gauss.getW(k) * mem.getDetJ(k);
+        qV -= fk * gradk * gauss.getW(k) * e.getDetJ(k);
     }
     return qV;
 }
@@ -138,7 +133,6 @@ Eigen::MatrixXd HeatTerm::computeMatrix(Element const &e, std::vector<double> co
     // Get precomputed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Gauss integration
     Eigen::MatrixXd out = Eigen::Matrix2d::Zero(); // 2 dimensions, should match fk size
@@ -148,7 +142,7 @@ Eigen::MatrixXd HeatTerm::computeMatrix(Element const &e, std::vector<double> co
         Eigen::MatrixXd fk;
         f.eval(e, u, k, fk, false);
 
-        out += fk * gauss.getW(k) * mem.getDetJ(k);
+        out += fk * gauss.getW(k) * e.getDetJ(k);
     }
     return out;
 }
diff --git a/heat/src/wProblem.cpp b/heat/src/wProblem.cpp
index 997214de6e308584081a5e694e89853d85906bba..717fd161feec3940a9cbc34e86ae9bade48f4913 100644
--- a/heat/src/wProblem.cpp
+++ b/heat/src/wProblem.cpp
@@ -20,7 +20,7 @@
 #include "wSource.h"
 #include "wPeriodic.h" // sinon:  warning C4150: deletion of pointer to incomplete type 'heat::Periodic'; no destructor called
 #include "wElement.h"
-#include "wMem.h"
+
 #include "wTag.h"
 using namespace heat;
 
@@ -39,21 +39,20 @@ Problem::~Problem()
 }
 
 /**
- * @brief Update the elements memory (Jacobian)
- * @authors Adrien Crovato
+ * @brief Initialize the elements precomputed values
  */
-void Problem::updateMem()
+void Problem::initElems()
 {
     // Update volume Jacobian
     for (auto vol : media)
         for (auto e : vol->tag->elems)
-            e->getVMem().update(true);
+            e->initValues(true);
     // Update volume Jacobian (volume source)
     for (auto vol : srcs)
         for (auto e : vol->tag->elems)
-            e->getVMem().update(true);
+            e->initValues(true);
     // Update surface Jacobian (Neumann B.C.)
     for (auto surf : bnds)
         for (auto e : surf->tag->elems)
-            e->getVMem().update(false);
+            e->initValues(false);
 }
\ No newline at end of file
diff --git a/heat/src/wProblem.h b/heat/src/wProblem.h
index dd385f9330aab5387810a358a2f0133a7bd3bf69..c961be28815f692a64cfe723a112fa59b8a34ff8 100644
--- a/heat/src/wProblem.h
+++ b/heat/src/wProblem.h
@@ -57,7 +57,7 @@ public:
     void add(std::shared_ptr<Periodic> p) { pdic = p; }
 
 #ifndef SWIG
-    void updateMem();
+    void initElems();
     virtual void write(std::ostream &out) const override;
 #endif
 };
diff --git a/heat/src/wSolver.cpp b/heat/src/wSolver.cpp
index 685cabd12719e02da28920275a06d8cc604eee26..35bab22c0acd6493f3769b0cfbcc82630ed702c8 100644
--- a/heat/src/wSolver.cpp
+++ b/heat/src/wSolver.cpp
@@ -62,7 +62,7 @@ Solver::Solver(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver
     save = true;
 
     // update element memory
-    pbl->updateMem();
+    pbl->initElems();
 }
 
 Solver::~Solver()
@@ -266,7 +266,7 @@ void Solver::start(MshExport *mshWriter)
                              // (flux moyen . volume) sur l'element
                              Eigen::VectorXd qV = HeatTerm::computeFlux(*e, T1, mat->k);
 
-                             double V = e->computeV();
+                             double V = e->getVol();
 
                              if (save)
                              {
diff --git a/katoptron/src/LinearSolver.h b/katoptron/src/LinearSolver.h
index 7cc27fc87b4ff5abeb774fcf441d389287c5c9d9..0fe02a2f7abee9daca881b2a74049d7a73ee7e48 100644
--- a/katoptron/src/LinearSolver.h
+++ b/katoptron/src/LinearSolver.h
@@ -8,7 +8,7 @@
 #include "wTag.h"
 #include "wMedium.h"
 #include "wNode.h"
-#include "wMem.h"
+
 #include "Teuchos_ParameterList.hpp"
 
 #include "wProblem.h"
diff --git a/katoptron/src/ResultsDecl.hpp b/katoptron/src/ResultsDecl.hpp
index 4eeef4aa176b68606844f56323d9ed82927fa2b4..9fac353203749ce64ad6867142c4912a433a5e81 100644
--- a/katoptron/src/ResultsDecl.hpp
+++ b/katoptron/src/ResultsDecl.hpp
@@ -8,7 +8,7 @@
 #include "wTag.h"
 #include "wMedium.h"
 #include "wNode.h"
-#include "wMem.h"
+
 #include "Teuchos_ParameterList.hpp"
 
 #include "wProblem.h"
diff --git a/katoptron/src/ResultsDef.hpp b/katoptron/src/ResultsDef.hpp
index b4e5ab91ad1994485e4f72c12210317822cb6fc6..985cfa1331676ba1fdd9a6b1e67a9f1494df6f4e 100644
--- a/katoptron/src/ResultsDef.hpp
+++ b/katoptron/src/ResultsDef.hpp
@@ -25,9 +25,9 @@ void writeResultsVTK(
     for (auto e : pbl.msh->elems)
     {
         if (e->type() == tbox::ELTYPE::HEX8 || e->type() == tbox::ELTYPE::TETRA4) // 3D
-            e->getVMem().update(true);
+            e->initValues(true);
         //else // 2D
-        //    e->getVMem().update(false);
+        //    e->initValues(false);
     }
     /* end temporary fix*/
 
diff --git a/katoptron/src/StressComputation.h b/katoptron/src/StressComputation.h
index 2870bd25e747f4c11792e58cf3a9b8140600d9d8..ec6c1e817215c04e83e718422b7c4a96f86b1d54 100644
--- a/katoptron/src/StressComputation.h
+++ b/katoptron/src/StressComputation.h
@@ -5,13 +5,11 @@
 #include "wSfHex8.h"
 #include "wGaussHex8.h"
 #include "wCacheHex8.h"
-#include "wMemHex8.h"
 
 #include "wTetra4.h"
 #include "wSfTetra4.h"
 #include "wGaussTetra4.h"
 #include "wCacheTetra4.h"
-#include "wMemTetra4.h"
 
 /**
  * @brief Those functions compute the stress and strain at nodes and in average per element.
@@ -19,9 +17,9 @@
  * such that it will be possible to use parallel for in the future.
  */
 
-void strain_at_gp(size_t n_n, tbox::Mem &mem, tbox::Cache &cache, size_t a, const Eigen::MatrixXd &u_e, Eigen::MatrixXd &Strain)
+void strain_at_gp(size_t n_n, tbox::Element &el, tbox::Cache &cache, size_t a, const Eigen::MatrixXd &u_e, Eigen::MatrixXd &Strain)
 {
-    Eigen::MatrixXd const &J = mem.getJinv(a);
+    Eigen::MatrixXd const &J = el.getJinv(a);
     Eigen::MatrixXd const &dffi = cache.getDsf(a);
     Eigen::Matrix3d du_d_Phi = Eigen::Matrix3d::Zero();
     for (size_t i = 0; i < n_n; ++i)
@@ -57,7 +55,6 @@ void compute_stress(size_t n_n, const Eigen::MatrixXd &Strain, Eigen::MatrixXd &
 void strain_stress_x_y_z(tbox::Element *elem, Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector<double> &x, const std::vector<double> &y, const std::vector<double> &z, const std::vector<double> &T, bool thermal = false)
 {
     tbox::Cache &cache = elem->getVCache();
-    tbox::Mem &mem = elem->getVMem();
 
     size_t n_gp = cache.getVGauss().getN();
     size_t n_n = elem->nodes.size();
@@ -74,7 +71,7 @@ void strain_stress_x_y_z(tbox::Element *elem, Eigen::MatrixXd &Strain, Eigen::Ma
     Strain = Eigen::Matrix3d::Zero();
     for (size_t a = 0; a < n_gp; ++a)
     {
-        strain_at_gp(n_n, mem, cache, a, u_e, Strain);
+        strain_at_gp(n_n, *elem, cache, a, u_e, Strain);
     }
     Strain /= n_n;
 
@@ -84,7 +81,6 @@ void strain_stress_x_y_z(tbox::Element *elem, Eigen::MatrixXd &Strain, Eigen::Ma
 void strain_stress_x_y_z_at_node(tbox::Element *elem, std::vector<Eigen::MatrixXd> &Strain_at_node, std::vector<Eigen::MatrixXd> &Stress_at_node, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector<double> &x, const std::vector<double> &y, const std::vector<double> &z, const std::vector<double> &T, bool thermal = false)
 {
     tbox::Cache &cache = elem->getVCache();
-    tbox::Mem &mem = elem->getVMem();
 
     size_t n_gp = cache.getVGauss().getN();
     size_t n_n = elem->nodes.size();
@@ -109,7 +105,7 @@ void strain_stress_x_y_z_at_node(tbox::Element *elem, std::vector<Eigen::MatrixX
     for (size_t a = 0; a < n_gp; ++a)
     {
         Strain = Eigen::Matrix3d::Zero();
-        strain_at_gp(n_n, mem, cache, a, u_e, Strain);
+        strain_at_gp(n_n, *elem, cache, a, u_e, Strain);
         Strain_at_gp[a] = Strain;
         compute_stress(n_n, Strain, Stress_at_gp[a], H, D, T, thermal);
     }
diff --git a/mirrors/src/wProblem.cpp b/mirrors/src/wProblem.cpp
index a708feeda80b0a3dcddc78a5f2bd4628485e8db8..16bd8f1363984404b49410423acef0b9c71b4b6d 100644
--- a/mirrors/src/wProblem.cpp
+++ b/mirrors/src/wProblem.cpp
@@ -23,7 +23,7 @@
 #include "wTNeumann.h"
 #include "wTSource.h"
 #include "wElement.h"
-#include "wMem.h"
+
 #include "wTag.h"
 using namespace mirrors;
 
@@ -32,26 +32,25 @@ Problem::Problem(std::shared_ptr<MshData> _msh, std::string _problem_type, doubl
 }
 
 /**
- * @brief Update the elements memory (Jacobian)
- * @authors Adrien Crovato
+ * @brief Initialize the elements precomputed values
  */
-void Problem::updateMem()
+void Problem::initElems()
 {
     // Update volume Jacobian
     for (auto vol : media)
         for (auto e : vol->tag->elems)
-            e->getVMem().update(true);
+            e->initValues(true);
     // Update volume Jacobian (volumic sources)
     for (auto vol : Tsrcs)
         for (auto e : vol->tag->elems)
-            e->getVMem().update(true);
+            e->initValues(true);
     // Update surface Jacobian (Neumann B.C.)
     for (auto surf : nuBCs)
         for (auto e : surf->tag->elems)
-            e->getVMem().update(false);
+            e->initValues(false);
     for (auto surf : nTBCs)
         for (auto e : surf->tag->elems)
-            e->getVMem().update(false);
+            e->initValues(false);
 }
 
 void Problem::write(std::ostream &out) const
diff --git a/mirrors/src/wProblem.h b/mirrors/src/wProblem.h
index 21a61be696d4e52d2a39e572a21186374c4e7eeb..0ce6de2b5911cc328a3e67ee4eda5b951a21b69c 100644
--- a/mirrors/src/wProblem.h
+++ b/mirrors/src/wProblem.h
@@ -62,7 +62,7 @@ public:
     Problem(std::shared_ptr<MshData> _msh, std::string problem_type, double T_ref = 0);
 
 #ifndef SWIG
-    void updateMem();
+    void initElems();
     virtual void write(std::ostream &out) const override;
 #endif
 };
diff --git a/mirrors/src/wSolver.cpp b/mirrors/src/wSolver.cpp
index ef21f944c25d3136485a14c69951b80b2d782b83..7a01d5c41dd90b38e7a206246a3cea1054cff0b5 100644
--- a/mirrors/src/wSolver.cpp
+++ b/mirrors/src/wSolver.cpp
@@ -50,7 +50,7 @@
 #include <sstream>
 #include <iostream>
 /*
-// [RB] not needed anymore 
+// [RB] not needed anymore
 #ifdef max
 #undef max    // pour std::numeric_limits<std::streamsize>::max()
 #endif
@@ -62,7 +62,7 @@ using namespace mirrors;
 Solver::Solver(Problem &_pbl, std::shared_ptr<tbox::LinearSolver> _linsol, int _nthreads) : pbl(_pbl), linsol(_linsol), nthreads(_nthreads)
 {
     // Update element memory
-    pbl.updateMem();
+    pbl.initElems();
 }
 
 void Solver::write(std::ostream &out) const
@@ -620,7 +620,7 @@ void Solver::start(MshExport *mshWriter)
 
             Tag &current_tag = *(ms->tag);
 
-            Eigen::Vector3d normal = current_tag.elems[0]->normal();
+            Eigen::Vector3d normal = current_tag.elems[0]->getNormal();
 
             double max_distance = 0.;
 
diff --git a/mirrors/src/wThermoMecaTerm.cpp b/mirrors/src/wThermoMecaTerm.cpp
index d8a005e514e73d069390e4e06f0a6c1803aa6447..712542b2d697d6db39ae819f124cc929a8289670 100644
--- a/mirrors/src/wThermoMecaTerm.cpp
+++ b/mirrors/src/wThermoMecaTerm.cpp
@@ -19,7 +19,6 @@
 #include "wElement.h"
 #include "wCache.h"
 #include "wGauss.h"
-#include "wMem.h"
 #include "wNode.h"
 
 using namespace mirrors;
@@ -33,7 +32,6 @@ Eigen::MatrixXd ThermoMecaTerm::buildM(Element const &e)
     // Get precomputed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Mass matrix
     Eigen::MatrixXd M = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
@@ -41,7 +39,7 @@ Eigen::MatrixXd ThermoMecaTerm::buildM(Element const &e)
     {
         // Shape functions
         Eigen::VectorXd const &ff = cache.getSf(k);
-        M += ff * ff.transpose() * gauss.getW(k) * mem.getDetJ(k);
+        M += ff * ff.transpose() * gauss.getW(k) * e.getDetJ(k);
     }
     return M;
 }
@@ -54,7 +52,6 @@ Eigen::MatrixXd ThermoMecaTerm::buildK(tbox::Element const &e, Eigen::MatrixXd c
     // Get precomputed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
     size_t ns = e.nodes.size(); // number of nodes
 
     // Meca matrix
@@ -62,7 +59,7 @@ Eigen::MatrixXd ThermoMecaTerm::buildK(tbox::Element const &e, Eigen::MatrixXd c
     for (size_t k = 0; k < gauss.getN(); ++k)
     {
         // Jacobian inverse
-        Eigen::Matrix3d const &invJ = mem.getJinv(k);
+        Eigen::Matrix3d const &invJ = e.getJinv(k);
         // Fill B matrix (shape functions)
         Eigen::MatrixXd const &dff = cache.getDsf(k);
         Eigen::MatrixXd B = Eigen::MatrixXd::Zero(6, 3 * ns);
@@ -82,7 +79,7 @@ Eigen::MatrixXd ThermoMecaTerm::buildK(tbox::Element const &e, Eigen::MatrixXd c
         }
 
         // Compute stiffness matrix
-        K += B.transpose() * H * B * gauss.getW(k) * mem.getDetJ(k);
+        K += B.transpose() * H * B * gauss.getW(k) * e.getDetJ(k);
     }
     return K;
 }
@@ -95,15 +92,14 @@ Eigen::MatrixXd ThermoMecaTerm::buildL(Element const &e, Eigen::MatrixXd const &
     // Get precomputed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Thermal matrix
     Eigen::MatrixXd L = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
     for (size_t k = 0; k < gauss.getN(); ++k)
     {
         // Jacobian
-        double detJ = mem.getDetJ(k);
-        Eigen::Matrix3d const &J = mem.getJinv(k);
+        double detJ = e.getDetJ(k);
+        Eigen::Matrix3d const &J = e.getJinv(k);
         Eigen::MatrixXd const &dff = cache.getDsf(k);
 
         // [AC] I refactored this from an entry by entry assembly to a full matrix multiplication,
@@ -121,7 +117,6 @@ Eigen::MatrixXd ThermoMecaTerm::buildS(Element const &e, Eigen::MatrixXd const &
     // Get precomputed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
     size_t ns = e.nodes.size(); // number of nodes
 
     // TM matrix
@@ -129,8 +124,8 @@ Eigen::MatrixXd ThermoMecaTerm::buildS(Element const &e, Eigen::MatrixXd const &
     for (size_t k = 0; k < gauss.getN(); ++k)
     {
         // Jacobian and shape functions
-        double detJ = mem.getDetJ(k);
-        Eigen::Matrix3d const &J = mem.getJinv(k);
+        double detJ = e.getDetJ(k);
+        Eigen::Matrix3d const &J = e.getJinv(k);
         Eigen::MatrixXd const &dff = cache.getDsf(k);
 
         // S
@@ -148,7 +143,6 @@ void ThermoMecaTerm::strain_stress(Element const &e, Eigen::MatrixXd const &H, s
     // Get precomputed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
     size_t ns = e.nodes.size(); // number of nodes
 
     Strain = Stress = Eigen::Matrix3d::Zero(); // 3 dimensions
@@ -163,7 +157,7 @@ void ThermoMecaTerm::strain_stress(Element const &e, Eigen::MatrixXd const &H, s
     for (size_t k = 0; k < gauss.getN(); ++k)
     {
         // Jacobian and strain
-        Eigen::Matrix3d const &J = mem.getJinv(k);
+        Eigen::Matrix3d const &J = e.getJinv(k);
         Eigen::MatrixXd const &dff = cache.getDsf(k);
         Eigen::Matrix3d du_d_Phi = Eigen::Matrix3d::Zero();
         for (size_t i = 0; i < ns; ++i)
diff --git a/tbox/_src/tboxw.h b/tbox/_src/tboxw.h
index 85f8528942b6cae7fb4440cbd522fe2df4d9b6f6..b9feacdafc52341b820365f648c374ef1b68ced9 100644
--- a/tbox/_src/tboxw.h
+++ b/tbox/_src/tboxw.h
@@ -43,7 +43,6 @@
 #include "wLazy.h"
 #include "wLine2.h"
 #include "wLinearSolver.h"
-#include "wMemTri3.h"
 #include "wMshConvert.h"
 #include "wMshCrack.h"
 #include "wMshData.h"
diff --git a/tbox/src/tbox.h b/tbox/src/tbox.h
index 4c2304c6db28aadfdd01ead77a77c8e36499c3be..b6814d4c6444cc2afffd02c09a043b438db500b7 100644
--- a/tbox/src/tbox.h
+++ b/tbox/src/tbox.h
@@ -44,7 +44,6 @@ class Tag;
 
 // Element handling
 class Cache;
-class Mem;
 
 // Point1
 class Point1;
@@ -54,35 +53,30 @@ class CachePoint1;
 class SfLine2;
 class CacheLine2;
 class GaussLine2;
-class MemLine2;
 class Line2;
 
 // Tri3
 class SfTri3;
 class CacheTri3;
 class GaussTri3;
-class MemTri3;
 class Tri3;
 
 // Quad4
 class SfQuad4;
 class CacheQuad4;
 class GaussQuad4;
-class MemQuad4;
 class Quad4;
 
 // Tetra4
 class SfTetra4;
 class CacheTetra4;
 class GaussTetra4;
-class MemTetra4;
 class Tetra4;
 
 // Hex8
 class SfHex8;
 class CacheHex8;
 class GaussHex8;
-class MemHex8;
 class Hex8;
 
 // Mesh utilities
diff --git a/tbox/src/wElement.cpp b/tbox/src/wElement.cpp
index aa8fba025e73710a0ced04b034a157f469686070..2aab2292503540e5df27b0a66c32fc4dd0dd9abe 100644
--- a/tbox/src/wElement.cpp
+++ b/tbox/src/wElement.cpp
@@ -17,8 +17,6 @@
 #include "wElement.h"
 #include "wNode.h"
 #include "wTag.h"
-#include "wFct0.h"
-#include "wMem.h"
 #include "wCache.h"
 #include "std_extra.h"
 
@@ -61,8 +59,8 @@ operator<<(std::ostream &out, ELTYPE const &obj)
 
 Element::Element(int n, Tag *_ptag, Tag *_etag,
                  std::vector<Tag *> &_parts, std::vector<Node *> &nods)
-    : wObject(), order(2), no(n), ptag(_ptag), etag(_etag), parts(_parts),
-      nodes(nods)
+    : wObject(), order(0), hasVals(false), hasGrads(false), isVol(false),
+      no(n), ptag(_ptag), etag(_etag), parts(_parts), nodes(nods)
 {
     if (ptag)
         ptag->elems.push_back(this);
@@ -73,37 +71,69 @@ Element::Element(int n, Tag *_ptag, Tag *_etag,
         p->elems.push_back(this);
 }
 
-Mem &Element::getVMem() const
+void Element::initValues(bool isvol, int ordr)
 {
-    throw std::runtime_error("Element::getVMem not implemented\n");
+    throw std::runtime_error("Element::initValues() not implemented!\n");
+}
+void Element::initGradients()
+{
+    throw std::runtime_error("Element::initGradients() not implemented!\n");
+}
+void Element::update()
+{
+    throw std::runtime_error("Element::update() not implemented!\n");
 }
 
-double Element::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
+/**
+ * @brief Check that determinant of Jacobian is positive
+ */
+void Element::checkJac(size_t k) const
 {
-    throw std::runtime_error("Element::buildJ not implemented\n");
+    if (detJs[k] <= 0)
+    {
+        std::stringstream err;
+        err << "Negative jacobian detected for element: " << *this << "\n";
+        throw std::runtime_error(err.str());
+    }
 }
-double Element::buildJ(size_t k) const
+
+void Element::buildJ()
 {
     throw std::runtime_error("Element::buildJ not implemented\n");
 }
-void Element::buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eigen::MatrixXd> &dJ) const
+void Element::buildSurfaceJ()
 {
-    throw std::runtime_error("Element::buildGradientJ not implemented\n");
+    throw std::runtime_error("Element::buildSurfaceJ not implemented\n");
 }
-void Element::buildGradientJ(size_t k, std::vector<double> &dDetJ) const
+void Element::buildGradientJ()
 {
     throw std::runtime_error("Element::buildGradientJ not implemented\n");
 }
+void Element::buildGradientSurfaceJ()
+{
+    throw std::runtime_error("Element::buildGradientSurfaceJ not implemented\n");
+}
 
-double Element::computeV() const
+void Element::computeV()
 {
     throw std::runtime_error("Element::computeV not implemented\n");
 }
-std::vector<double> Element::buildGradientV(int dim) const
+void Element::buildGradientV(int dim)
 {
     throw std::runtime_error("Element::buildGradientV not implemented\n");
 }
 
+/**
+ * @brief Compute center of gravity of element
+ */
+void Element::computeCg()
+{
+    cg = Eigen::Vector3d::Zero();
+    for (auto n : nodes)
+        cg += n->pos;
+    cg /= nodes.size();
+}
+
 std::vector<Eigen::Vector3d> Element::computeTangents() const
 {
     throw std::runtime_error("Element::computeTangents not implemented\n");
@@ -113,13 +143,13 @@ std::vector<Eigen::Vector3d> Element::buildTangents(size_t k) const
     throw std::runtime_error("Element::buildTangents not implemented\n");
 }
 
-Eigen::Vector3d Element::normal() const
+void Element::computeNormal()
 {
-    throw std::runtime_error("Element::normal not implemented\n");
+    throw std::runtime_error("Element::computeNormal not implemented\n");
 }
-std::vector<Eigen::Vector3d> Element::gradientNormal() const
+void Element::computeGradientNormal()
 {
-    throw std::runtime_error("Element::gradientNormal not implemented\n");
+    throw std::runtime_error("Element::computeGradientNormal not implemented\n");
 }
 
 /**
diff --git a/tbox/src/wElement.h b/tbox/src/wElement.h
index 72c4dd9b92680ff801c2583ab7c8fa8ef41eeccc..3c02267a31885b6d3f8bcc961503387628244f38 100644
--- a/tbox/src/wElement.h
+++ b/tbox/src/wElement.h
@@ -49,8 +49,27 @@ TBOX_API std::ostream &operator<<(std::ostream &out, ELTYPE const &obj);
  */
 class TBOX_API Element : public fwk::wObject
 {
+    friend class MshCrack; // so that MshCrack can access normal and cg before update
 protected:
     int order; ///< order of integration (Gauss quadrature)
+    // Flags
+    bool hasVals;  ///< true if (precomputed) values have been initialized
+    bool hasGrads; ///< true if (precomputed) gradients have been initialized
+    bool isVol;    ///< true if element is a volume element
+    // Jacobian
+    std::vector<double> detJs;                        ///< determinant of Jacobian matrix
+    std::vector<Eigen::MatrixXd> Js;                  ///< Jacobian matrix
+    std::vector<Eigen::MatrixXd> iJs;                 ///< inverse of Jacobian matrix
+    std::vector<std::vector<double>> gradDetJs;       ///< gradient of Jacobian determinant wrt to nodes position
+    std::vector<std::vector<Eigen::MatrixXd>> gradJs; ///< gradient of Jacobian matrix wrt to nodes position
+    // Volume and CG
+    double vol;                  ///< element volume
+    Eigen::Vector3d cg;          ///< element center of gravity
+    std::vector<double> gradVol; ///< gradient of volume wrt to nodes position
+    // Normal vectors
+    Eigen::Vector3d normal;                  ///< unit normal vector
+    std::vector<Eigen::Vector3d> gradNormal; ///< graident of unit normal vector wrt to nodes position
+
 public:
     int no;                    ///< label
     Tag *ptag;                 ///< physical tag   (group)
@@ -70,30 +89,49 @@ public:
 #ifndef SWIG
 protected:
     // Jacobian
-    virtual double buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const;
-    virtual double buildJ(size_t k) const;
-    virtual void buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eigen::MatrixXd> &dJ) const;
-    virtual void buildGradientJ(size_t k, std::vector<double> &dDetJ) const;
+    void checkJac(size_t k) const;
+    virtual void buildJ();
+    virtual void buildSurfaceJ();
+    virtual void buildGradientJ();
+    virtual void buildGradientSurfaceJ();
+    // CG and Volume
+    void computeCg();
+    virtual void computeV();
+    virtual void buildGradientV(int dim);
+    // Tangent and normal vectors
+    virtual std::vector<Eigen::Vector3d> computeTangents() const;
+    virtual std::vector<Eigen::Vector3d> buildTangents(size_t k) const;
+    virtual void computeNormal();
+    virtual void computeGradientNormal();
 
 public:
+    // Initialize and update precomputed values
+    virtual void initValues(bool isvol, int ordr = 2);
+    virtual void initGradients();
+    virtual void update();
     // Getters
     virtual Cache &getVCache() const = 0;
-    virtual Mem &getVMem() const;
-    // Volume (@todo private)
-    virtual double computeV() const;
-    virtual std::vector<double> buildGradientV(int dim) const;
-    // Tangent and normal vectors (@todo private)
-    virtual std::vector<Eigen::Vector3d> computeTangents() const;
-    virtual std::vector<Eigen::Vector3d> buildTangents(size_t k) const;
-    virtual Eigen::Vector3d normal() const;
-    virtual std::vector<Eigen::Vector3d> gradientNormal() const;
+    inline bool hasValues() const;
+    inline bool hasGradients() const;
+    inline double getDetJ(size_t k) const;
+    inline Eigen::MatrixXd const &getJ(size_t k) const;
+    inline Eigen::MatrixXd const &getJinv(size_t k) const;
+    inline std::vector<double> const &getGradDetJ(size_t k) const;
+    inline std::vector<Eigen::MatrixXd> const &getGradJ(size_t k) const;
+    inline double getVol() const;
+    inline Eigen::Vector3d const &getCg() const;
+    inline std::vector<double> const &getGradVol() const;
+    inline Eigen::Vector3d const &getNormal() const;
+    inline std::vector<Eigen::Vector3d> const &getGradNormal() const;
     // Utilities
     virtual Eigen::VectorXd computeGradient(std::vector<double> const &u, size_t k) const;
-    virtual Eigen::MatrixXd gp2node(Eigen::MatrixXd &values_at_gp) const;
+    Eigen::MatrixXd gp2node(Eigen::MatrixXd &values_at_gp) const;
     virtual void write(std::ostream &out) const override;
 #endif
 };
 
+#include "wElement.inl"
+
 } // namespace tbox
 
 #endif //WELEMENT_H
diff --git a/tbox/src/wElement.inl b/tbox/src/wElement.inl
new file mode 100644
index 0000000000000000000000000000000000000000..9bed1207dc092d4da6b4cbc58f601dd91421dc0e
--- /dev/null
+++ b/tbox/src/wElement.inl
@@ -0,0 +1,166 @@
+/*
+ * Copyright 2020 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/**
+ * @brief Return wether the element values have been initialized or not
+ */
+inline bool Element::hasValues() const
+{
+    return hasVals;
+}
+
+/**
+ * @brief Return wether the element gradients have been initialized or not
+ */
+inline bool Element::hasGradients() const
+{
+    return hasGrads;
+}
+
+/**
+ * @brief Return the determinant of the Jacobian at Gauss point k
+ */
+inline double Element::getDetJ(size_t k) const
+{
+    if (hasVals)
+        return detJs[k];
+    else
+        throw std::runtime_error("Element::getDetJ is invalid because element memory is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the the Jacobian matrix at Gauss point k
+ */
+inline Eigen::MatrixXd const &Element::getJ(size_t k) const
+{
+    if (hasVals)
+    {
+        if (isVol)
+            return Js[k];
+        else
+            throw std::runtime_error("Element::getJ Jacobian matrix is not available for surface element!\n");
+    }
+    else
+        throw std::runtime_error("Element::getJ is invalid because element memory is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the inverse of the Jacobian matrix at Gauss point k
+ */
+inline Eigen::MatrixXd const &Element::getJinv(size_t k) const
+{
+    if (hasVals)
+    {
+        if (isVol)
+            return iJs[k];
+        else
+            throw std::runtime_error("Element::getJinv inverse Jacobian matrix is not available for surface element!\n");
+    }
+    else
+        throw std::runtime_error("Element::getJinv is invalid because element memory is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the volume of the element
+ */
+inline double Element::getVol() const
+{
+    if (hasVals)
+        return vol;
+    else
+        throw std::runtime_error("Mem::getVol is invalid because element memory is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the center of gravity of the element
+ */
+inline Eigen::Vector3d const &Element::getCg() const
+{
+    if (hasVals)
+        return cg;
+    else
+        throw std::runtime_error("Element::getCg is invalid because element memory is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the unit normal vector of the element
+ */
+inline Eigen::Vector3d const &Element::getNormal() const
+{
+    if (hasVals)
+    {
+        if (!isVol)
+            return normal;
+        else
+            throw std::runtime_error("Element::getNormal normal vector is not available for volume element!\n");
+    }
+    else
+        throw std::runtime_error("Element::getNormal is invalid because element memory is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the gradient of determinant of the Jacobian with respect to node positions at Gauss point k
+ */
+inline std::vector<double> const &Element::getGradDetJ(size_t k) const
+{
+    if (hasGrads)
+        return gradDetJs[k];
+    else
+        throw std::runtime_error("Element::getGradDetJ is invalid because element memory (gradient) is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the gradient of the Jacobian matrix with respect to node positions at Gauss point k
+ */
+inline std::vector<Eigen::MatrixXd> const &Element::getGradJ(size_t k) const
+{
+    if (hasGrads)
+    {
+        if (isVol)
+            return gradJs[k];
+        else
+            throw std::runtime_error("Element::getGradJ gradient of Jacobian matrix is not available for surface element!\n");
+    }
+    else
+        throw std::runtime_error("Element::getGradJ is invalid because element memory (gradient) is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the gradient of the volume with respect to node positions
+ */
+inline std::vector<double> const &Element::getGradVol() const
+{
+    if (hasGrads)
+        return gradVol;
+    else
+        throw std::runtime_error("Element::getGradVol is invalid because element memory (gradient) is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the gradient of the unit normal vector with respect to node positions
+ */
+inline std::vector<Eigen::Vector3d> const &Element::getGradNormal() const
+{
+    if (hasGrads)
+    {
+        if (!isVol)
+            return gradNormal;
+        else
+            throw std::runtime_error("Element::getGradNormal normal vector gradient is not available for volume element!\n");
+    }
+    else
+        throw std::runtime_error("Element::getGradNormal is invalid because element memory (gradient) is not up-to-date!\n");
+}
diff --git a/tbox/src/wHex8.cpp b/tbox/src/wHex8.cpp
index ae4c523af2c6a7ad95f3fa13b3667711ec455eac..606ce1063b5857afc6a1437182152181cf4070d7 100644
--- a/tbox/src/wHex8.cpp
+++ b/tbox/src/wHex8.cpp
@@ -18,7 +18,6 @@
 #include "wSfHex8.h"
 #include "wGaussHex8.h"
 #include "wCacheHex8.h"
-#include "wMemHex8.h"
 #include "wNode.h"
 
 using namespace tbox;
@@ -26,60 +25,108 @@ using namespace tbox;
 Hex8::Hex8(int n, Tag *_ptag, Tag *_etag,
            std::vector<Tag *> &_parts, std::vector<Node *> &nods) : Element(n, _ptag, _etag, _parts, nods)
 {
-    pmem = new MemHex8(*this, static_cast<int>(getCache(order).gauss.getN()));
 }
-Hex8::~Hex8()
+
+/**
+ * @brief Return cache handling shape functions and Gauss points for all Hex8 elements
+ */
+Cache &Hex8::getVCache() const
 {
-    if (pmem)
-        delete pmem;
+    return getCache(order);
 }
 
 /**
- * @brief Return private memory handling Jacobian for this Hex8 element 
+ * @brief Initialize precomputed values
  */
-Mem &Hex8::getVMem() const
+void Hex8::initValues(bool isvol, int ordr)
 {
-    return getMem();
+    // Set flags
+    isVol = isvol;
+    order = ordr;
+    hasVals = true;
+    // Update values
+    size_t ngp = getCache(order).getVGauss().getN();
+    Js.resize(ngp);
+    iJs.resize(ngp);
+    detJs.resize(ngp);
+    this->update();
 }
-
 /**
- * @brief Return cache handling shape functions and Gauss points for all Hex8 elements
+ * @brief Initialize precomputed gradients
  */
-Cache &Hex8::getVCache() const
+void Hex8::initGradients()
 {
-    return getCache(order);
+    // Set flags
+    hasGrads = true;
+    // Update values
+    size_t ngp = getCache(order).getVGauss().getN();
+    gradJs.resize(ngp);
+    gradDetJs.resize(ngp);
+    this->update();
+}
+/**
+ * @brief Update Hex8 precomputed values and gradients
+ */
+void Hex8::update()
+{
+    // Update values
+    if (hasVals)
+    {
+        // Jacobian
+        if (!isVol)
+            throw std::runtime_error("Hex8::update(): the problem must always have the same dimension as Hex8!\n");
+        // Hex8 is always of same dimension as Problem (volume element)
+        this->buildJ();
+        // Volume and CG
+        this->computeV();
+        this->computeCg();
+    }
+    // Update gradients
+    if (hasGrads)
+    {
+        // Jacobian
+        if (!isVol)
+            throw std::runtime_error("Hex8::updateGradients(): the problem must always have the same dimension as Hex8!\n");
+        // Hex8 is always of same dimension as Problem (volume element)
+        this->buildGradientJ();
+        // Volume
+        this->buildGradientV(3);
+    }
 }
 
 /**
- * @brief Compute Jacobian matrix, inverse and determinant at Gauss point k
- * 
+ * @brief Compute Jacobian matrix, inverse and determinant
+ *
  * J_ij(3,3) = diN_n * xj_n (n is the node)
  */
-double Hex8::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
+void Hex8::buildJ()
 {
-    Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
+    CacheHex8 &cache = getCache(order);
+    for (size_t k = 0; k < cache.getVGauss().getN(); ++k)
+    {
+        Eigen::MatrixXd const &dff = cache.getDsf(k);
 
-    Eigen::Matrix3d JJ = Eigen::Matrix3d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
-    size_t i = 0;
-    for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
-        JJ += dff.col(i) * (*it)->pos.transpose();
-    J = JJ;
-    iJ = JJ.inverse();
-    return JJ.determinant();
+        Eigen::Matrix3d JJ = Eigen::Matrix3d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
+        size_t i = 0;
+        for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
+            JJ += dff.col(i) * (*it)->pos.transpose();
+        Js[k] = JJ;
+        iJs[k] = JJ.inverse();
+        detJs[k] = JJ.determinant();
+        this->checkJac(k);
+    }
 }
 
 /**
  * @brief Compute the element volume
  */
-double Hex8::computeV() const
+void Hex8::computeV()
 {
     CacheHex8 &cache = getCache(order);
     GaussHex8 &gauss = cache.gauss;
-    MemHex8 &mem = getMem();
 
     // Gauss integration
-    double V = 0.0;
+    vol = 0.0;
     for (size_t k = 0; k < gauss.getN(); ++k)
-        V += gauss.getW(k) * mem.getDetJ(k);
-    return V;
+        vol += gauss.getW(k) * this->getDetJ(k);
 }
diff --git a/tbox/src/wHex8.h b/tbox/src/wHex8.h
index bf39aaa5123a926d59a7182cb9d8da7b329b0a32..391b73dce4b71ae7b5866f7d2b26dfc8bf1fe19b 100644
--- a/tbox/src/wHex8.h
+++ b/tbox/src/wHex8.h
@@ -32,28 +32,24 @@ namespace tbox
 class TBOX_API Hex8 : public Element
 {
 #ifndef SWIG
-    friend class MemHex8;
-
 private:
-    MemHex8 *pmem; ///< private memory of precalculated values
     inline static CacheHex8 &getCache(int n);
-    inline MemHex8 &getMem() const;
 
 protected:
-    virtual double buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const override;
-
-public: // @todo
-    virtual double computeV() const override;
+    virtual void buildJ() override;
+    virtual void computeV() override;
 #endif
 
 public:
     Hex8(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
-    virtual ~Hex8();
+    virtual ~Hex8() {}
     virtual ELTYPE type() const override { return ELTYPE::HEX8; }
 
 #ifndef SWIG
+    virtual void initValues(bool isvol, int ordr = 2) override;
+    virtual void initGradients() override;
+    virtual void update() override;
     virtual Cache &getVCache() const override;
-    virtual Mem &getVMem() const override;
 #endif
 };
 
diff --git a/tbox/src/wHex8.inl b/tbox/src/wHex8.inl
index 6912e483f986152c6952c5bc7ef2c45372ebe2e5..a2c43adf6081419b973c51cf173693b702155969 100644
--- a/tbox/src/wHex8.inl
+++ b/tbox/src/wHex8.inl
@@ -14,14 +14,6 @@
  * limitations under the License.
  */
 
-/**
- * @brief Return private memory for Hex8 private methods
- */
-inline MemHex8 &Hex8::getMem() const
-{
-    return *pmem;
-}
-
 /**
  * @brief Return cache for Hex8 private methods
  */
diff --git a/tbox/src/wLine2.cpp b/tbox/src/wLine2.cpp
index 85dce33b026da3146b4f99bda70d84fd5bd29a1d..b28d131676469869e3835a7a57ea1b40bfcdfe31 100644
--- a/tbox/src/wLine2.cpp
+++ b/tbox/src/wLine2.cpp
@@ -18,7 +18,6 @@
 #include "wSfLine2.h"
 #include "wGaussLine2.h"
 #include "wCacheLine2.h"
-#include "wMemLine2.h"
 #include "wNode.h"
 
 using namespace tbox;
@@ -26,61 +25,114 @@ using namespace tbox;
 Line2::Line2(int n, Tag *_ptag, Tag *_etag,
              std::vector<Tag *> &_parts, std::vector<Node *> &nods) : Element(n, _ptag, _etag, _parts, nods)
 {
-    pmem = new MemLine2(*this, static_cast<int>(getCache(order).gauss.getN()));
 }
-Line2::~Line2()
+
+/**
+ * @brief Return cache holding shape functions and Gauss points for all Line2 elements
+ */
+Cache &Line2::getVCache() const
 {
-    if (pmem)
-        delete pmem;
+    return getCache(order);
 }
 
 /**
- * @brief Return private memory handling Jacobian for this Line2 element
+ * @brief Initialize precomputed values
  */
-Mem &Line2::getVMem() const
+void Line2::initValues(bool isvol, int ordr)
 {
-    return getMem();
+    // Set flags
+    isVol = isvol;
+    order = ordr;
+    hasVals = true;
+    // Update values
+    size_t ngp = getCache(order).getVGauss().getN();
+    Js.resize(ngp);
+    iJs.resize(ngp);
+    detJs.resize(ngp);
+    this->update();
 }
-
 /**
- * @brief Return cache holding shape functions and Gauss points for all Line2 elements
+ * @brief Initialize precomputed gradients
  */
-Cache &Line2::getVCache() const
+void Line2::initGradients()
 {
-    return getCache(order);
+    // Set flags
+    hasGrads = true;
+    // Update values
+    size_t ngp = getCache(order).getVGauss().getN();
+    gradJs.resize(ngp);
+    gradDetJs.resize(ngp);
+    this->update();
+}
+/**
+ * @brief Update precomputed values and gradients
+ */
+void Line2::update()
+{
+    // Update values
+    if (hasVals)
+    {
+        // Jacobian
+        if (isVol)
+            throw std::runtime_error("Line2::update(): the problem must always have a dimension less than Line2!\n");
+        // Line2 is never of same dimension as Problem (surface element)
+        this->buildSurfaceJ();
+        // Length and CG
+        this->computeV();
+        this->computeCg();
+        // Normal
+        this->computeNormal();
+    }
+    // Update gradients
+    if (hasGrads)
+    {
+        // Jacobian
+        if (isVol)
+            throw std::runtime_error("Line2::updateGradient(): the problem must always have a dimension less than Line2!\n");
+        // Line2 is never of same dimension as Problem (surface element)
+        this->buildGradientSurfaceJ();
+        // Length
+        this->buildGradientV(2);
+        // Normal
+        this->computeGradientNormal();
+    }
 }
 
 /**
- * @brief Compute Jacobian determinant at Gauss point k
+ * @brief Compute Jacobian determinant
  *
  * detJ = norm(x_i * dN_i)
  */
-double Line2::buildJ(size_t k) const
+void Line2::buildSurfaceJ()
 {
-    return this->buildTangents(k)[0].norm();
+    for (size_t k = 0; k < getCache(order).getVGauss().getN(); ++k)
+        detJs[k] = this->buildTangents(k)[0].norm();
 }
 
 /**
- * @brief Compute Jacobian determinant gradient with respect to each node coordinate, at Gauss point k
+ * @brief Compute Jacobian determinant gradient with respect to each node coordinate
  *
  * dDetJ[l = i * 2 + j] = detJ_l * dt0_l;
  */
-void Line2::buildGradientJ(size_t k, std::vector<double> &dDetJ) const
+void Line2::buildGradientSurfaceJ()
 {
-    Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
-
-    // Tangent vectors and determinant
-    Eigen::Vector3d detJ = this->buildTangents(k)[0].normalized();
-
-    // Jacobian gradient
-    dDetJ.resize(2 * 2); // 2 nodes and 2 dimensions (x, y)
-    for (size_t i = 0; i < 2; ++i)
+    CacheLine2 &cache = getCache(order);
+    for (size_t k = 0; k < cache.getVGauss().getN(); ++k)
     {
-        for (size_t j = 0; j < 2; ++j)
+        // Shape functions, tangent vectors and determinant
+        Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
+        Eigen::Vector3d detJ = this->buildTangents(k)[0].normalized();
+
+        // Jacobian gradient
+        gradDetJs[k].resize(2 * 2); // 2 nodes and 2 dimensions (x, y)
+        for (size_t i = 0; i < 2; ++i)
         {
-            Eigen::Vector3d dt0 = Eigen::Vector3d::Zero();
-            dt0(j) = dsf(0, i); // retain sf gradients associated to ith node for dimension j
-            dDetJ[i * 2 + j] = detJ.dot(dt0);
+            for (size_t j = 0; j < 2; ++j)
+            {
+                Eigen::Vector3d dt0 = Eigen::Vector3d::Zero();
+                dt0(j) = dsf(0, i); // retain sf gradients associated to ith node for dimension j
+                gradDetJs[k][i * 2 + j] = detJ.dot(dt0);
+            }
         }
     }
 }
@@ -88,38 +140,34 @@ void Line2::buildGradientJ(size_t k, std::vector<double> &dDetJ) const
 /**
  * @brief Compute the element length
  */
-double Line2::computeV() const
+void Line2::computeV()
 {
     CacheLine2 &cache = getCache(order);
     GaussLine2 &gauss = cache.gauss;
-    MemLine2 &mem = getMem();
 
     // Gauss integration
-    double V = 0.0;
+    vol = 0.0;
     for (size_t k = 0; k < gauss.getN(); ++k)
-        V += gauss.getW(k) * mem.getDetJ(k);
-    return V;
+        vol += gauss.getW(k) * this->getDetJ(k);
 }
 /**
  * @brief Compute length gradient with respect to each node coordinate
  *
  * dV[i] = w_k * dDetJ[i]_k
  */
-std::vector<double> Line2::buildGradientV(int dim) const
+void Line2::buildGradientV(int dim)
 {
     GaussLine2 &gauss = getCache(order).gauss;
-    MemLine2 &mem = getMem();
 
     // Gradient of volume
-    std::vector<double> dVol(2 * dim, 0.); // 2 nodes
+    gradVol.resize(2 * dim, 0.); // 2 nodes
     for (size_t k = 0; k < gauss.getN(); ++k)
     {
-        double w = gauss.getW(k);                       // GP weight
-        std::vector<double> dDetJ = mem.getGradDetJ(k); // gradient of Jacobian determinant
+        double w = gauss.getW(k);                         // GP weight
+        std::vector<double> dDetJ = this->getGradDetJ(k); // gradient of Jacobian determinant
         for (size_t i = 0; i < dDetJ.size(); ++i)
-            dVol[i] += w * dDetJ[i];
+            gradVol[i] += w * dDetJ[i];
     }
-    return dVol;
 }
 
 /**
@@ -149,16 +197,16 @@ std::vector<Eigen::Vector3d> Line2::buildTangents(size_t k) const
 /**
  * @brief Compute element unit normal vector
  */
-Eigen::Vector3d Line2::normal() const
+void Line2::computeNormal()
 {
     Eigen::Vector3d t = this->computeTangents()[0].normalized();
-    return Eigen::Vector3d(-t(1), t(0), 0.);
+    normal = Eigen::Vector3d(-t(1), t(0), 0.);
 }
 
 /**
  * @brief Compute gradient of element unit normal vector with respect to each node coordinate
  */
-std::vector<Eigen::Vector3d> Line2::gradientNormal() const
+void Line2::computeGradientNormal()
 {
     // Compute normal and norm
     Eigen::Vector3d t = this->computeTangents()[0];
@@ -167,10 +215,9 @@ std::vector<Eigen::Vector3d> Line2::gradientNormal() const
 
     // Compute gradients
     Eigen::Matrix3d inn = (Eigen::Matrix3d::Identity() - n * n.transpose() / (nn * nn)) / nn; // gradient of n/|n|
-    std::vector<Eigen::Vector3d> dn(2 * 2);                                                   // 2 nodes and 2 dimensions
-    dn[0] = inn * Eigen::Vector3d(0, -1, 0);
-    dn[1] = inn * Eigen::Vector3d(1, 0, 0);
-    dn[2] = inn * Eigen::Vector3d(0, 1, 0);
-    dn[3] = inn * Eigen::Vector3d(-1, 0, 0); // outer gradient times gradient of t (pos) wtr x
-    return dn;
+    gradNormal.resize(2 * 2);                                                                 // 2 nodes and 2 dimensions
+    gradNormal[0] = inn * Eigen::Vector3d(0, -1, 0);
+    gradNormal[1] = inn * Eigen::Vector3d(1, 0, 0);
+    gradNormal[2] = inn * Eigen::Vector3d(0, 1, 0);
+    gradNormal[3] = inn * Eigen::Vector3d(-1, 0, 0); // outer gradient times gradient of t (pos) wtr x
 }
diff --git a/tbox/src/wLine2.h b/tbox/src/wLine2.h
index 9e6e5f0de9dc05c7a676de610640515aec8b1baa..6b076608d104fc2ffe5877ee05ed9a69c918d68a 100644
--- a/tbox/src/wLine2.h
+++ b/tbox/src/wLine2.h
@@ -32,33 +32,29 @@ namespace tbox
 class TBOX_API Line2 : public Element
 {
 #ifndef SWIG
-    friend class MemLine2;
-
 private:
-    MemLine2 *pmem; ///< private memory of precalculated values
-    inline MemLine2 &getMem() const;
     inline static CacheLine2 &getCache(int n);
 
 protected:
-    virtual double buildJ(size_t k) const override;
-    virtual void buildGradientJ(size_t k, std::vector<double> &dDetJ) const override;
-
-public: // @todo
-    virtual double computeV() const override;
-    virtual std::vector<double> buildGradientV(int dim) const override;
-    virtual Eigen::Vector3d normal() const override;
-    virtual std::vector<Eigen::Vector3d> gradientNormal() const override;
+    virtual void buildSurfaceJ() override;
+    virtual void buildGradientSurfaceJ() override;
+    virtual void computeV() override;
+    virtual void buildGradientV(int dim) override;
+    virtual void computeNormal() override;
+    virtual void computeGradientNormal() override;
     virtual std::vector<Eigen::Vector3d> computeTangents() const override;
     virtual std::vector<Eigen::Vector3d> buildTangents(size_t k) const override;
 #endif
 
 public:
     Line2(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
-    virtual ~Line2();
+    virtual ~Line2() {}
     virtual ELTYPE type() const override { return ELTYPE::LINE2; }
 
 #ifndef SWIG
-    virtual Mem &getVMem() const override;
+    virtual void initValues(bool isvol, int ordr = 2) override;
+    virtual void initGradients() override;
+    virtual void update() override;
     virtual Cache &getVCache() const override;
 #endif
 };
diff --git a/tbox/src/wLine2.inl b/tbox/src/wLine2.inl
index 704f23146100a2e1ac7115e688f82e3ee85f6205..5dc091848a30ebf69769361d298fe45e83b04408 100644
--- a/tbox/src/wLine2.inl
+++ b/tbox/src/wLine2.inl
@@ -14,14 +14,6 @@
  * limitations under the License.
  */
 
-/**
- * @brief Return private memory for Line2 private methods
- */
-inline MemLine2 &Line2::getMem() const
-{
-    return *pmem;
-}
-
 /**
  * @brief Return cache for Line2 private methods
  */
diff --git a/tbox/src/wMem.cpp b/tbox/src/wMem.cpp
deleted file mode 100644
index 4a003a0c6edd42121ee83548021b95c66d9eceaf..0000000000000000000000000000000000000000
--- a/tbox/src/wMem.cpp
+++ /dev/null
@@ -1,59 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wMem.h"
-#include "wElement.h"
-
-using namespace tbox;
-
-Mem::Mem(size_t _ngp) : valid(false), validG(false), ngp(_ngp), detJs(_ngp), Js(_ngp), Jinvs(_ngp)
-{
-    vol = 0;
-    cg = Eigen::Vector3d::Zero();
-}
-
-/**
- * @brief Check that determinant of Jacobian is positive
- */
-void Mem::checkJac(Element const &e, double detJ) const
-{
-    if (detJ <= 0)
-    {
-        std::stringstream err;
-        err << "Negative jacobian detected for element: " << e << "\n";
-        throw std::runtime_error(err.str());
-    }
-}
-
-/**
- * @brief Initialize the gradients
- */
-void Mem::initGradients(bool sameDim)
-{
-    gradJs.resize(ngp);
-    gradDetJs.resize(ngp);
-    this->updateGradients(sameDim);
-}
-
-void Mem::update(bool sameDim)
-{
-    throw std::runtime_error("Mem::update() not implemented!\n");
-}
-
-void Mem::updateGradients(bool sameDim)
-{
-    throw std::runtime_error("Mem::updateGradients() not implemented!\n");
-}
diff --git a/tbox/src/wMem.h b/tbox/src/wMem.h
deleted file mode 100644
index 8a8945f4479f475a63c03ff640028dd5c8448aed..0000000000000000000000000000000000000000
--- a/tbox/src/wMem.h
+++ /dev/null
@@ -1,86 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WMEM_H
-#define WMEM_H
-
-#include "tbox.h"
-#include <vector>
-#include <Eigen/Dense>
-
-#ifndef SWIG
-
-namespace tbox
-{
-
-/**
- * @brief Base class for element private memory
- * @authors Romain Boman, Adrien Crovato
- * @todo Remove capability to getJ and getJinv when not volume element
- * @todo Split values and gradients in 2 different Mem?
- */
-class TBOX_API Mem
-{
-protected:
-    // Flag
-    bool valid;  ///< validity flag (true if memory has been updated once)
-    bool validG; ///< validity flag for gradients (true if gradients have been initialized)
-
-    // Element value
-    size_t ngp;         ///< number of Gauss points
-    double vol;         ///< element volume
-    Eigen::Vector3d cg; ///< element center of gravity
-
-    // GP values
-    std::vector<double> detJs;          ///< determinant of Jacobian matrix
-    std::vector<Eigen::MatrixXd> Js;    ///< Jacobian matrix
-    std::vector<Eigen::MatrixXd> Jinvs; ///< inverse of Jacobian matrix
-
-    // Gradients
-    std::vector<double> gradVol;                      ///< gradient of volume wrt to nodes position
-    std::vector<std::vector<double>> gradDetJs;       ///< gradient of Jacobian determinant wrt to nodes position
-    std::vector<std::vector<Eigen::MatrixXd>> gradJs; ///< gradient of Jacobian matrix wrt to nodes position
-
-    // Functions
-    void checkJac(Element const &e, double detJ) const;
-    virtual void updateGradients(bool sameDim);
-
-public:
-    explicit Mem(size_t _ngp);
-    virtual ~Mem() {}
-
-    inline bool isValid();
-    inline bool isGradientValid();
-    inline double getDetJ(size_t k);
-    inline Eigen::MatrixXd const &getJ(size_t k);
-    inline Eigen::MatrixXd const &getJinv(size_t k);
-    inline double getVol();
-    inline Eigen::Vector3d const &getCg();
-    inline std::vector<double> const &getGradDetJ(size_t k);
-    inline std::vector<Eigen::MatrixXd> const &getGradJ(size_t k);
-    inline std::vector<double> const &getGradVol();
-
-    void initGradients(bool sameDim);
-    virtual void update(bool sameDim);
-};
-
-#include "wMem.inl"
-
-} // namespace tbox
-
-#endif
-
-#endif //WMEM_H
diff --git a/tbox/src/wMem.inl b/tbox/src/wMem.inl
deleted file mode 100644
index 217af82ab85bb06859d1e84cf863141e94d4580b..0000000000000000000000000000000000000000
--- a/tbox/src/wMem.inl
+++ /dev/null
@@ -1,119 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-/**
- * @brief Return the validity flag
- */
-inline bool Mem::isValid()
-{
-    return valid;
-}
-
-/**
- * @brief Return the validity flag (gradients)
- */
-inline bool Mem::isGradientValid()
-{
-    return validG;
-}
-
-/**
- * @brief Return the determinant of the Jacobian at Gauss point k
- */
-inline double Mem::getDetJ(size_t k)
-{
-    if (valid)
-        return detJs[k];
-    else
-        throw std::runtime_error("Mem::getDetJ is invalid because element memory is not up-to-date!\n");
-}
-
-/**
- * @brief Return the the Jacobian matrix at Gauss point k
- */
-inline Eigen::MatrixXd const &Mem::getJ(size_t k)
-{
-    if (valid)
-        return Js[k];
-    else
-        throw std::runtime_error("Mem::getJ is invalid because element memory is not up-to-date!\n");
-}
-
-/**
- * @brief Return the inverse of the Jacobian matrix at Gauss point k
- */
-inline Eigen::MatrixXd const &Mem::getJinv(size_t k)
-{
-    if (valid)
-        return Jinvs[k];
-    else
-        throw std::runtime_error("Mem::getJinv is invalid because element memory is not up-to-date!\n");
-}
-
-/**
- * @brief Return the volume of the element
- */
-inline double Mem::getVol()
-{
-    if (valid)
-        return vol;
-    else
-        throw std::runtime_error("Mem::getVol is invalid because element memory is not up-to-date!\n");
-}
-
-/**
- * @brief Return the center of gravity of the element
- */
-inline Eigen::Vector3d const &Mem::getCg()
-{
-    if (valid)
-        return cg;
-    else
-        throw std::runtime_error("Mem::getCg is invalid because element memory is not up-to-date!\n");
-}
-
-/**
- * @brief Return the gradient of determinant of the Jacobian with respect to node positions at Gauss point k
- */
-inline std::vector<double> const &Mem::getGradDetJ(size_t k)
-{
-    if (validG)
-        return gradDetJs[k];
-    else
-        throw std::runtime_error("Mem::getGradDetJ is invalid because element memory (gradient) is not up-to-date!\n");
-}
-
-/**
- * @brief Return the gradient of the Jacobian matrix with respect to node positions at Gauss point k
- */
-inline std::vector<Eigen::MatrixXd> const &Mem::getGradJ(size_t k)
-{
-    if (validG)
-        return gradJs[k];
-    else
-        throw std::runtime_error("Mem::getGradJ is invalid because element memory (gradient) is not up-to-date!\n");
-}
-
-/**
- * @brief Return the gradient of the volume with respect to node positions
- */
-inline std::vector<double> const &Mem::getGradVol()
-{
-    if (validG)
-        return gradVol;
-    else
-        throw std::runtime_error("Mem::getGradVol is invalid because element memory (gradient) is not up-to-date!\n");
-}
diff --git a/tbox/src/wMemHex8.cpp b/tbox/src/wMemHex8.cpp
deleted file mode 100644
index c815f2025241d25ffaf9c35f40d079999f1c0fdd..0000000000000000000000000000000000000000
--- a/tbox/src/wMemHex8.cpp
+++ /dev/null
@@ -1,68 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wMemHex8.h"
-#include "wNode.h"
-
-using namespace tbox;
-
-MemHex8::MemHex8(Hex8 &_e, size_t _ngp) : Mem(_ngp), e(_e)
-{
-}
-
-/**
- * @brief Update Hex8 values
- */
-void MemHex8::update(bool sameDim)
-{
-    // Jacobian
-    if (!sameDim)
-        throw std::runtime_error("MemHex8::update(): the problem must always have the same dimension as Hex8!\n");
-    // Hex8 is always of same dimension as Problem (volume element)
-    for (int k = 0; k < ngp; ++k)
-    {
-        detJs[k] = e.buildJ(k, Js[k], Jinvs[k]);
-        checkJac(e, detJs[k]);
-    }
-    // CG
-    cg = Eigen::Vector3d::Zero();
-    for (auto n : e.nodes)
-        cg += n->pos;
-    cg /= e.nodes.size();
-    // Memory is now valid and the volume can be computed
-    valid = true;
-    vol = e.computeV();
-
-    // Update gradient if they are needed
-    if (validG)
-        this->updateGradients(sameDim);
-}
-
-/**
- * @brief Update Hex8 gradients
- */
-void MemHex8::updateGradients(bool sameDim)
-{
-    // Jacobian
-    if (!sameDim)
-        throw std::runtime_error("MemHex8::updateGradients(): the problem must always have the same dimension as Hex8!\n");
-    // Hex8 is always of same dimension as Problem (volume element)
-    for (int k = 0; k < ngp; ++k)
-        e.buildGradientJ(k, gradDetJs[k], gradJs[k]);
-    // Volume
-    validG = true;
-    gradVol = e.buildGradientV(3);
-}
diff --git a/tbox/src/wMemHex8.h b/tbox/src/wMemHex8.h
deleted file mode 100644
index 0787c7bbdb9aea926e845701f7e24cdd9c18f560..0000000000000000000000000000000000000000
--- a/tbox/src/wMemHex8.h
+++ /dev/null
@@ -1,50 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WMEMHEX8_H
-#define WMEMHEX8_H
-
-#include "tbox.h"
-#include "wMem.h"
-#include "wHex8.h"
-
-#ifndef SWIG
-
-namespace tbox
-{
-
-/**
- * @brief Precalculated values for Hex8
- * @authors Romain Boman, Adrien Crovato
- */
-class TBOX_API MemHex8 : public Mem
-{
-    // Element type
-    Hex8 &e;
-
-public:
-    explicit MemHex8(Hex8 &_e, size_t _ngp);
-    virtual ~MemHex8() {}
-
-    virtual void update(bool sameDim) override;
-    virtual void updateGradients(bool sameDim) override;
-};
-
-} // namespace tbox
-
-#endif
-
-#endif //WMEMHEX8_H
diff --git a/tbox/src/wMemLine2.cpp b/tbox/src/wMemLine2.cpp
deleted file mode 100644
index 8e0c8022ec509f26c9e55bb14802cbb62368c011..0000000000000000000000000000000000000000
--- a/tbox/src/wMemLine2.cpp
+++ /dev/null
@@ -1,65 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wMemLine2.h"
-#include "wNode.h"
-
-using namespace tbox;
-
-MemLine2::MemLine2(Line2 &_e, size_t _ngp) : Mem(_ngp), e(_e)
-{
-}
-
-/**
- * @brief Update Line2 values
- */
-void MemLine2::update(bool sameDim)
-{
-    // Jacobian
-    if (sameDim)
-        throw std::runtime_error("MemLine2::update(): the problem must always have a dimension less than Line2!\n");
-    // Line2 is never of same dimension as Problem (surface element)
-    for (int k = 0; k < ngp; ++k)
-        detJs[k] = e.buildJ(k);
-    // CG
-    cg = Eigen::Vector3d::Zero();
-    for (auto n : e.nodes)
-        cg += n->pos;
-    cg /= e.nodes.size();
-    // Memory is now valid and the volume can be computed
-    valid = true;
-    vol = e.computeV();
-
-    // Update gradient if they are needed
-    if (validG)
-        this->updateGradients(sameDim);
-}
-
-/**
- * @brief Update Line2 gradients
- */
-void MemLine2::updateGradients(bool sameDim)
-{
-    // Jacobian
-    if (sameDim)
-        throw std::runtime_error("MemLine2::updateGradient(): the problem must always have a dimension less than Line2!\n");
-    // Line2 is never of same dimension as Problem (surface element)
-    for (int k = 0; k < ngp; ++k)
-        e.buildGradientJ(k, gradDetJs[k]);
-    // Volume
-    validG = true;
-    gradVol = e.buildGradientV(2);
-}
diff --git a/tbox/src/wMemLine2.h b/tbox/src/wMemLine2.h
deleted file mode 100644
index 38624af08e2a5d6a9ac369ae7487ea55ae47aa10..0000000000000000000000000000000000000000
--- a/tbox/src/wMemLine2.h
+++ /dev/null
@@ -1,50 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WMEMLINE2_H
-#define WMEMLINE2_H
-
-#include "tbox.h"
-#include "wMem.h"
-#include "wLine2.h"
-
-#ifndef SWIG
-
-namespace tbox
-{
-
-/**
- * @brief Precalculated values for Line2
- * @authors Romain Boman, Adrien Crovato
- */
-class TBOX_API MemLine2 : public Mem
-{
-    // Element type
-    Line2 &e;
-
-public:
-    explicit MemLine2(Line2 &_e, size_t _ngp);
-    virtual ~MemLine2() {}
-
-    virtual void update(bool sameDim) override;
-    virtual void updateGradients(bool sameDim) override;
-};
-
-} // namespace tbox
-
-#endif
-
-#endif //WMEMLINE2_H
diff --git a/tbox/src/wMemQuad4.cpp b/tbox/src/wMemQuad4.cpp
deleted file mode 100644
index 50ee3b51eeb5e3fceed1e16ac27ccf2f6e4b2b63..0000000000000000000000000000000000000000
--- a/tbox/src/wMemQuad4.cpp
+++ /dev/null
@@ -1,81 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wMemQuad4.h"
-#include "wNode.h"
-
-using namespace tbox;
-
-MemQuad4::MemQuad4(Quad4 &_e, size_t _ngp) : Mem(_ngp), e(_e)
-{
-}
-
-/**
- * @brief Update Quad4 values
- */
-void MemQuad4::update(bool sameDim)
-{
-    // Jacobian
-    // Quad4 is not of same dimension as Problem (surface element)
-    if (!sameDim)
-    {
-        for (int k = 0; k < ngp; ++k)
-            detJs[k] = e.buildJ(k);
-    }
-    // Quad4 is of same dimension as Problem (volume element)
-    else
-    {
-        for (int k = 0; k < ngp; ++k)
-        {
-            detJs[k] = e.buildJ(k, Js[k], Jinvs[k]);
-            checkJac(e, detJs[k]);
-        }
-    }
-    // CG
-    cg = Eigen::Vector3d::Zero();
-    for (auto n : e.nodes)
-        cg += n->pos;
-    cg /= e.nodes.size();
-    // Memory is now valid and the volume can be computed
-    valid = true;
-    vol = e.computeV();
-
-    // Update gradient if they are needed
-    if (validG)
-        this->updateGradients(sameDim);
-}
-
-/**
- * @brief Update Quad4 gradients
- */
-void MemQuad4::updateGradients(bool sameDim)
-{
-    // Jacobian
-    // Quad4 is not of same dimension as Problem (surface element)
-    if (!sameDim)
-        for (int k = 0; k < ngp; ++k)
-            e.buildGradientJ(k, gradDetJs[k]);
-    // Quad4 is of same dimension as Problem (volume element)
-    else
-        for (int k = 0; k < ngp; ++k)
-            e.buildGradientJ(k, gradDetJs[k], gradJs[k]);
-    // Volume
-    validG = true;
-    if (!sameDim)
-        gradVol = e.buildGradientV(3); // gradient of surface
-    else
-        gradVol = e.buildGradientV(2); // gradient of volume
-}
diff --git a/tbox/src/wMemQuad4.h b/tbox/src/wMemQuad4.h
deleted file mode 100644
index cfe653f1c6dbfec444764cd17d5e9c219ad4d682..0000000000000000000000000000000000000000
--- a/tbox/src/wMemQuad4.h
+++ /dev/null
@@ -1,50 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WMEMQUAD4_H
-#define WMEMQUAD4_H
-
-#include "tbox.h"
-#include "wMem.h"
-#include "wQuad4.h"
-
-#ifndef SWIG
-
-namespace tbox
-{
-
-/**
- * @brief Precalculated values for Quad4
- * @authors Romain Boman, Adrien Crovato
- */
-class TBOX_API MemQuad4 : public Mem
-{
-    // Element type
-    Quad4 &e;
-
-public:
-    explicit MemQuad4(Quad4 &_e, size_t _ngp);
-    virtual ~MemQuad4() {}
-
-    virtual void update(bool sameDim) override;
-    virtual void updateGradients(bool sameDim) override;
-};
-
-} // namespace tbox
-
-#endif
-
-#endif //WMEMQUAD4_H
diff --git a/tbox/src/wMemTetra4.cpp b/tbox/src/wMemTetra4.cpp
deleted file mode 100644
index 5be98b0c5ab2aeccd9dd2b28ad58899c27588df2..0000000000000000000000000000000000000000
--- a/tbox/src/wMemTetra4.cpp
+++ /dev/null
@@ -1,68 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wMemTetra4.h"
-#include "wNode.h"
-
-using namespace tbox;
-
-MemTetra4::MemTetra4(Tetra4 &_e, size_t _ngp) : Mem(_ngp), e(_e)
-{
-}
-
-/**
- * @brief Update Tetra4 values
- */
-void MemTetra4::update(bool sameDim)
-{
-    // Jacobian
-    if (!sameDim)
-        throw std::runtime_error("MemTetra4::update(): the problem must always have the same dimension as Tetra4!\n");
-    // Tetra4 is always of same dimension as Problem (volume element)
-    for (size_t k = 0; k < ngp; ++k)
-    {
-        detJs[k] = e.buildJ(k, Js[k], Jinvs[k]);
-        checkJac(e, detJs[k]);
-    }
-    // CG
-    cg = Eigen::Vector3d::Zero();
-    for (auto n : e.nodes)
-        cg += n->pos;
-    cg /= e.nodes.size();
-    // Memory is now valid and the volume can be computed
-    valid = true;
-    vol = e.computeV();
-
-    // Update gradient if they are needed
-    if (validG)
-        this->updateGradients(sameDim);
-}
-
-/**
- * @brief Update Tetra4 gradients
- */
-void MemTetra4::updateGradients(bool sameDim)
-{
-    // Jacobian
-    if (!sameDim)
-        throw std::runtime_error("MemTetra4::updateGradients(): the problem must always have the same dimension as Tetra4!\n");
-    // Tetra4 is always of same dimension as Problem (volume element)
-    for (int k = 0; k < ngp; ++k)
-        e.buildGradientJ(k, gradDetJs[k], gradJs[k]);
-    // Volume
-    validG = true;
-    gradVol = e.buildGradientV(3);
-}
diff --git a/tbox/src/wMemTetra4.h b/tbox/src/wMemTetra4.h
deleted file mode 100644
index 42a71dcbd34f63f7090d92674e4e16e7558092a4..0000000000000000000000000000000000000000
--- a/tbox/src/wMemTetra4.h
+++ /dev/null
@@ -1,50 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WMEMTETRA4_H
-#define WMEMTETRA4_H
-
-#include "tbox.h"
-#include "wMem.h"
-#include "wTetra4.h"
-
-#ifndef SWIG
-
-namespace tbox
-{
-
-/**
- * @brief Precalculated values for Tetra4
- * @authors Romain Boman, Adrien Crovato
- */
-class TBOX_API MemTetra4 : public Mem
-{
-    // Element type
-    Tetra4 &e;
-
-public:
-    explicit MemTetra4(Tetra4 &_e, size_t _ngp);
-    virtual ~MemTetra4() {}
-
-    virtual void update(bool sameDim) override;
-    virtual void updateGradients(bool sameDim) override;
-};
-
-} // namespace tbox
-
-#endif
-
-#endif //WMEMTETRA4_H
diff --git a/tbox/src/wMemTri3.cpp b/tbox/src/wMemTri3.cpp
deleted file mode 100644
index 707e7d55838171ad620a1c20f85f5f1b7547aaaf..0000000000000000000000000000000000000000
--- a/tbox/src/wMemTri3.cpp
+++ /dev/null
@@ -1,81 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wMemTri3.h"
-#include "wNode.h"
-
-using namespace tbox;
-
-MemTri3::MemTri3(Tri3 &_e, size_t _ngp) : Mem(_ngp), e(_e)
-{
-}
-
-/**
- * @brief Update Tri3 values
- */
-void MemTri3::update(bool sameDim)
-{
-    // Jacobian
-    // Tri3 is not of same dimension as Problem (surface element)
-    if (!sameDim)
-    {
-        for (int k = 0; k < ngp; ++k)
-            detJs[k] = e.buildJ(k);
-    }
-    // Tri3 is of same dimension as Problem (volume element)
-    else
-    {
-        for (int k = 0; k < ngp; ++k)
-        {
-            detJs[k] = e.buildJ(k, Js[k], Jinvs[k]);
-            checkJac(e, detJs[k]);
-        }
-    }
-    // CG
-    cg = Eigen::Vector3d::Zero();
-    for (auto n : e.nodes)
-        cg += n->pos;
-    cg /= e.nodes.size();
-    // Memory is now valid and the volume can be computed
-    valid = true;
-    vol = e.computeV();
-
-    // Update gradient if they are needed
-    if (validG)
-        this->updateGradients(sameDim);
-}
-
-/**
- * @brief Update Tri3 gradients
- */
-void MemTri3::updateGradients(bool sameDim)
-{
-    // Jacobian
-    // Tri3 is not of same dimension as Problem (surface element)
-    if (!sameDim)
-        for (int k = 0; k < ngp; ++k)
-            e.buildGradientJ(k, gradDetJs[k]);
-    // Tri3 is of same dimension as Problem (volume element)
-    else
-        for (int k = 0; k < ngp; ++k)
-            e.buildGradientJ(k, gradDetJs[k], gradJs[k]);
-    // Volume
-    validG = true;
-    if (!sameDim)
-        gradVol = e.buildGradientV(3); // gradient of surface
-    else
-        gradVol = e.buildGradientV(2); // gradient of volume
-}
diff --git a/tbox/src/wMemTri3.h b/tbox/src/wMemTri3.h
deleted file mode 100644
index 025e54ea12f19ccce4abdbcf8d002e9f18d8476a..0000000000000000000000000000000000000000
--- a/tbox/src/wMemTri3.h
+++ /dev/null
@@ -1,50 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *     http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WMEMTRI3_H
-#define WMEMTRI3_H
-
-#include "tbox.h"
-#include "wMem.h"
-#include "wTri3.h"
-
-#ifndef SWIG
-
-namespace tbox
-{
-
-/**
- * @brief Precalculated values for Tri3
- * @authors Romain Boman, Adrien Crovato
- */
-class TBOX_API MemTri3 : public Mem
-{
-    // Element type
-    Tri3 &e;
-
-public:
-    explicit MemTri3(Tri3 &_e, size_t _ngp);
-    virtual ~MemTri3() {}
-
-    virtual void update(bool sameDim) override;
-    virtual void updateGradients(bool sameDim) override;
-};
-
-} // namespace tbox
-
-#endif
-
-#endif //WMEMTRI3_H
diff --git a/tbox/src/wMshCrack.cpp b/tbox/src/wMshCrack.cpp
index a54dc692697fb6375b912f140ca691586c8d1c00..8b550193a1099b2e1c12e8058a8bfaa3cb168571 100644
--- a/tbox/src/wMshCrack.cpp
+++ b/tbox/src/wMshCrack.cpp
@@ -225,11 +225,10 @@ void MshCrack::modifyBoundaries()
     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)
     {
-        crkCG[i] = Eigen::Vector3d::Zero();
-        for (size_t j = 0; j < crckGrp->tag->elems[i]->nodes.size(); ++j)
-            crkCG[i] += crckGrp->tag->elems[i]->nodes[j]->pos;
-        crkCG[i] /= crckGrp->tag->elems[i]->nodes.size();
-        crkN[i] = crckGrp->tag->elems[i]->normal();
+        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)
diff --git a/tbox/src/wMshDeform.cpp b/tbox/src/wMshDeform.cpp
index 6eacfebaa80367dea2cc0be5c5976513f3141a54..259fa27977259698875b23318ed3d683010ba81f 100644
--- a/tbox/src/wMshDeform.cpp
+++ b/tbox/src/wMshDeform.cpp
@@ -23,7 +23,6 @@
 #include "wElement.h"
 #include "wCache.h"
 #include "wGauss.h"
-#include "wMem.h"
 #include "wGmres.h"
 
 #include <deque>
@@ -241,7 +240,6 @@ Eigen::MatrixXd MshDeform::buildK(tbox::Element const &e, Eigen::MatrixXd const
     // Get precomputed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
     size_t ns = e.nodes.size(); // number of nodes
 
     // Elementary stiffness matrix
@@ -249,7 +247,7 @@ Eigen::MatrixXd MshDeform::buildK(tbox::Element const &e, Eigen::MatrixXd const
     for (size_t k = 0; k < gauss.getN(); ++k)
     {
         // Jacobian inverse
-        Eigen::MatrixXd const &iJ = mem.getJinv(k);
+        Eigen::MatrixXd const &iJ = e.getJinv(k);
         // Fill B matrix (shape functions)
         Eigen::MatrixXd const &dff = cache.getDsf(k);
         Eigen::MatrixXd B;
@@ -286,7 +284,7 @@ Eigen::MatrixXd MshDeform::buildK(tbox::Element const &e, Eigen::MatrixXd const
         }
 
         // Compute stiffness matrix
-        K += B.transpose() * H * B * gauss.getW(k) * mem.getDetJ(k);
+        K += B.transpose() * H * B * gauss.getW(k) * e.getDetJ(k);
     }
     return K;
 }
@@ -314,8 +312,8 @@ void MshDeform::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &K)
         if (nDim == 2)
         {
             H = Eigen::Matrix3d::Zero();
-            double nu = 0.;               // Poisson ratio [-1, 0.5]
-            double E = 1 / e->computeV(); // Young modulus
+            double nu = 0.;             // Poisson ratio [-1, 0.5]
+            double E = 1 / e->getVol(); // Young modulus
             H(0, 0) = 1 - nu;
             H(1, 1) = 1 - nu;
             H(2, 2) = (1 - 2 * nu) / 2;
@@ -326,8 +324,8 @@ void MshDeform::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &K)
         else
         {
             H = Eigen::MatrixXd::Zero(6, 6);
-            double nu = 0.;               // Poisson ratio [-1, 0.5]
-            double E = 1 / e->computeV(); // Young modulus
+            double nu = 0.;             // Poisson ratio [-1, 0.5]
+            double E = 1 / e->getVol(); // Young modulus
             H(0, 0) = 1 - nu;
             H(1, 1) = 1 - nu;
             H(2, 2) = 1 - nu;
@@ -494,13 +492,13 @@ void MshDeform::deform()
     }
     // Update element memory
     for (auto e : fldElems)
-        e->getVMem().update(true);
+        e->update();
     for (auto e : movBndElems)
-        e->getVMem().update(false);
+        e->update();
     for (auto e : symBndElems)
-        e->getVMem().update(false);
+        e->update();
     for (auto e : intBndElems)
-        e->getVMem().update(false);
+        e->update();
     std::cout << "done" << std::endl;
 }
 
diff --git a/tbox/src/wQuad4.cpp b/tbox/src/wQuad4.cpp
index d3981e8b5325e0dfd4b6b50447a65b2de5ae17fa..c2ee84f612d8a2343d12e7d2461dc3ef2f114cdb 100644
--- a/tbox/src/wQuad4.cpp
+++ b/tbox/src/wQuad4.cpp
@@ -18,7 +18,6 @@
 #include "wSfQuad4.h"
 #include "wGaussQuad4.h"
 #include "wCacheQuad4.h"
-#include "wMemQuad4.h"
 #include "wNode.h"
 
 using namespace tbox;
@@ -26,72 +25,135 @@ using namespace tbox;
 Quad4::Quad4(int n, Tag *_ptag, Tag *_etag,
              std::vector<Tag *> &_parts, std::vector<Node *> &nods) : Element(n, _ptag, _etag, _parts, nods)
 {
-    pmem = new MemQuad4(*this, static_cast<int>(getCache(order).gauss.getN()));
 }
-Quad4::~Quad4()
+
+/**
+ * @brief Return cache holding shape functions and Gauss points for all Quad4 elements
+ */
+Cache &Quad4::getVCache() const
 {
-    if (pmem)
-        delete pmem;
+    return getCache(order);
 }
 
 /**
- * @brief Return private memory handling Jacobian for this Quad4 element
+ * @brief Initialize precomputed values
  */
-Mem &Quad4::getVMem() const
+void Quad4::initValues(bool isvol, int ordr)
 {
-    return getMem();
+    // Set flags
+    isVol = isvol;
+    order = ordr;
+    hasVals = true;
+    // Update values
+    size_t ngp = getCache(order).getVGauss().getN();
+    Js.resize(ngp);
+    iJs.resize(ngp);
+    detJs.resize(ngp);
+    this->update();
 }
-
 /**
- * @brief Return cache holding shape functions and Gauss points for all Quad4 elements
+ * @brief Initialize precomputed gradients
  */
-Cache &Quad4::getVCache() const
+void Quad4::initGradients()
 {
-    return getCache(order);
+    // Set flags
+    hasGrads = true;
+    // Update values
+    size_t ngp = getCache(order).getVGauss().getN();
+    gradJs.resize(ngp);
+    gradDetJs.resize(ngp);
+    this->update();
+}
+/**
+ * @brief Update precomputed values and gradients
+ */
+void Quad4::update()
+{
+    // Update values
+    if (hasVals)
+    {
+        // Jacobian and normal
+        // Quad4 is not of same dimension as Problem (surface element)
+        if (!isVol)
+        {
+            this->buildSurfaceJ();
+            this->computeNormal();
+        }
+        // Quad4 is of same dimension as Problem (volume element)
+        else
+            this->buildJ();
+        // Volume and CG
+        this->computeV();
+        this->computeCg();
+    }
+    // Update gradients
+    if (hasGrads)
+    {
+        // Jacobian, volume and normal
+        // Quad4 is not of same dimension as Problem (surface element)
+        if (!isVol)
+        {
+            this->buildGradientSurfaceJ();
+            this->buildGradientV(3); // gradient of surface
+            this->computeGradientNormal();
+        }
+        // Quad4 is of same dimension as Problem (volume element)
+        else
+        {
+            this->buildGradientJ();
+            this->buildGradientV(2); // gradient of volume
+        }
+    }
 }
 
 /**
- * @brief Compute Jacobian determinant at Gauss point k
+ * @brief Compute Jacobian determinant
  *
  * detJ = norm(x_i dN_i,xi X x_i dN_i,eta)
  */
-double Quad4::buildJ(size_t k) const
+void Quad4::buildSurfaceJ()
 {
-    std::vector<Eigen::Vector3d> ts = this->buildTangents(k);
-    return ts[0].cross(ts[1]).norm();
+    for (size_t k = 0; k < getCache(order).getVGauss().getN(); ++k)
+    {
+        std::vector<Eigen::Vector3d> ts = this->buildTangents(k);
+        detJs[k] = ts[0].cross(ts[1]).norm();
+    }
 }
 /**
  * @brief Compute Jacobian matrix, inverse and determinant at Gauss point k
  *
  * J_ij(2,2) = diN_n * xj_n (n is the node)
  */
-double Quad4::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
+void Quad4::buildJ()
 {
-    Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
+    CacheQuad4 &cache = getCache(order);
+    for (size_t k = 0; k < cache.getVGauss().getN(); ++k)
+    {
+        Eigen::MatrixXd const &dff = cache.getDsf(k);
 
-    Eigen::Matrix2d JJ = Eigen::Matrix2d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
-    size_t i = 0;
-    for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
-        JJ += dff.col(i) * (*it)->pos.transpose().block<1, 2>(0, 0);
-    J = JJ;
-    iJ = JJ.inverse();
-    return JJ.determinant();
+        Eigen::Matrix2d JJ = Eigen::Matrix2d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
+        size_t i = 0;
+        for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
+            JJ += dff.col(i) * (*it)->pos.transpose().block<1, 2>(0, 0);
+        Js[k] = JJ;
+        iJs[k] = JJ.inverse();
+        detJs[k] = JJ.determinant();
+        this->checkJac(k);
+    }
 }
 
 /**
  * @brief Return the element surface
  */
-double Quad4::computeV() const
+void Quad4::computeV()
 {
     CacheQuad4 &cache = getCache(order);
     GaussQuad4 &gauss = cache.gauss;
-    MemQuad4 &mem = getMem();
 
     // Gauss integration
-    double V = 0.0;
+    vol = 0.0;
     for (size_t k = 0; k < gauss.getN(); ++k)
-        V += gauss.getW(k) * mem.getDetJ(k);
-    return V;
+        vol += gauss.getW(k) * this->getDetJ(k);
 }
 
 /**
@@ -126,8 +188,8 @@ std::vector<Eigen::Vector3d> Quad4::buildTangents(size_t k) const
 /**
  * @brief Compute element unit normal vector
  */
-Eigen::Vector3d Quad4::normal() const
+void Quad4::computeNormal()
 {
     std::vector<Eigen::Vector3d> ts = this->computeTangents();
-    return ts[0].cross(ts[1]).normalized();
+    normal = ts[0].cross(ts[1]).normalized();
 }
diff --git a/tbox/src/wQuad4.h b/tbox/src/wQuad4.h
index 0e3975b50c0d2323a7bb2f1c2bde585bc747c375..75cde412625a5cca2d1fb921bf209eef3e942006 100644
--- a/tbox/src/wQuad4.h
+++ b/tbox/src/wQuad4.h
@@ -32,32 +32,28 @@ namespace tbox
 class TBOX_API Quad4 : public Element
 {
 #ifndef SWIG
-    friend class MemQuad4;
-
 private:
-    MemQuad4 *pmem; ///< private memory of precalculated values
     inline static CacheQuad4 &getCache(int n);
-    inline MemQuad4 &getMem() const;
 
 protected:
-    virtual double buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const override;
-    virtual double buildJ(size_t k) const override;
-
-public: // @todo
-    virtual double computeV() const override;
+    virtual void buildJ() override;
+    virtual void buildSurfaceJ() override;
+    virtual void computeV() override;
     virtual std::vector<Eigen::Vector3d> computeTangents() const override;
     virtual std::vector<Eigen::Vector3d> buildTangents(size_t k) const override;
-    virtual Eigen::Vector3d normal() const override;
+    virtual void computeNormal() override;
 #endif
 
 public:
     Quad4(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
-    virtual ~Quad4();
+    virtual ~Quad4() {}
     virtual ELTYPE type() const override { return ELTYPE::QUAD4; }
 
 #ifndef SWIG
+    virtual void initValues(bool isvol, int ordr = 2) override;
+    virtual void initGradients() override;
+    virtual void update() override;
     virtual Cache &getVCache() const override;
-    virtual Mem &getVMem() const override;
 #endif
 };
 
diff --git a/tbox/src/wQuad4.inl b/tbox/src/wQuad4.inl
index 6e23c1b10e1aadfaa77f2b3663c933037199b74f..e77859ffbb6b5e75059c67da231cf2337604e0bd 100644
--- a/tbox/src/wQuad4.inl
+++ b/tbox/src/wQuad4.inl
@@ -14,14 +14,6 @@
  * limitations under the License.
  */
 
-/**
- * @brief Return private memory for Quad4 private methods
- */
-inline MemQuad4 &Quad4::getMem() const
-{
-    return *pmem;
-}
-
 /**
  * @brief Return cache for Quad4 private methods
  */
diff --git a/tbox/src/wTetra4.cpp b/tbox/src/wTetra4.cpp
index e7b8e5e3716cb6209013558bcc63c23cbd2b460a..8c1b42ab5a83974b2abe9e9e859b01eecb776ae5 100644
--- a/tbox/src/wTetra4.cpp
+++ b/tbox/src/wTetra4.cpp
@@ -18,7 +18,6 @@
 #include "wSfTetra4.h"
 #include "wGaussTetra4.h"
 #include "wCacheTetra4.h"
-#include "wMemTetra4.h"
 #include "wNode.h"
 
 using namespace tbox;
@@ -26,111 +25,160 @@ using namespace tbox;
 Tetra4::Tetra4(int n, Tag *_ptag, Tag *_etag,
                std::vector<Tag *> &_parts, std::vector<Node *> &nods) : Element(n, _ptag, _etag, _parts, nods)
 {
-    pmem = new MemTetra4(*this, static_cast<int>(getCache(order).gauss.getN()));
 }
-Tetra4::~Tetra4()
+
+/**
+ * @brief Return cache handling shape functions and Gauss points for all Tetra4 elements
+ */
+Cache &Tetra4::getVCache() const
 {
-    if (pmem)
-        delete pmem;
+    return getCache(order);
 }
 
 /**
- * @brief Return private memory handling Jacobian for this Tetra4 element
+ * @brief Initialize precomputed values
  */
-Mem &Tetra4::getVMem() const
+void Tetra4::initValues(bool isvol, int ordr)
 {
-    return getMem();
+    // Set flags
+    isVol = isvol;
+    order = ordr;
+    hasVals = true;
+    // Update values
+    size_t ngp = getCache(order).getVGauss().getN();
+    Js.resize(ngp);
+    iJs.resize(ngp);
+    detJs.resize(ngp);
+    this->update();
 }
-
 /**
- * @brief Return cache handling shape functions and Gauss points for all Tetra4 elements
+ * @brief Initialize precomputed gradients
  */
-Cache &Tetra4::getVCache() const
+void Tetra4::initGradients()
 {
-    return getCache(order);
+    // Set flags
+    hasGrads = true;
+    // Update values
+    size_t ngp = getCache(order).getVGauss().getN();
+    gradJs.resize(ngp);
+    gradDetJs.resize(ngp);
+    this->update();
+}
+/**
+ * @brief Update precomputed values and gradients
+ */
+void Tetra4::update()
+{
+    // Update gradients
+    if (hasVals)
+    {
+        // Jacobian
+        if (!isVol)
+            throw std::runtime_error("Tetra4::update(): the problem must always have the same dimension as Tetra4!\n");
+        // Tetra4 is always of same dimension as Problem (volume element)
+        this->buildJ();
+        // Volume and CG
+        this->computeV();
+        this->computeCg();
+    }
+    // Update gradients
+    if (hasGrads)
+    {
+        // Jacobian
+        if (!isVol)
+            throw std::runtime_error("Tetra4::updateGradients(): the problem must always have the same dimension as Tetra4!\n");
+        // Tetra4 is always of same dimension as Problem (volume element)
+        this->buildGradientJ();
+        // Volume
+        this->buildGradientV(3);
+    }
 }
 
 /**
- * @brief Compute Jacobian matrix, inverse and determinant at Gauss point k
+ * @brief Compute Jacobian matrix, inverse and determinant
  *
  * J_ij(3,3) = diN_n * xj_n (n is the node)
  */
-double Tetra4::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
+void Tetra4::buildJ()
 {
-    Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
+    CacheTetra4 &cache = getCache(order);
+    for (size_t k = 0; k < cache.getVGauss().getN(); ++k)
+    {
+        Eigen::MatrixXd const &dff = cache.getDsf(k);
 
-    Eigen::Matrix3d JJ = Eigen::Matrix3d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
-    size_t i = 0;
-    for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
-        JJ += dff.col(i) * ((*it)->pos).transpose();
-    J = JJ;
-    iJ = JJ.inverse();
-    return JJ.determinant();
+        Eigen::Matrix3d JJ = Eigen::Matrix3d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
+        size_t i = 0;
+        for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
+            JJ += dff.col(i) * ((*it)->pos).transpose();
+        Js[k] = JJ;
+        iJs[k] = JJ.inverse();
+        detJs[k] = JJ.determinant();
+        this->checkJac(k);
+    }
 }
 
 /**
- * @brief Compute Jacobian matrix gradient with respect to each node coordinate, at Gauss point k
+ * @brief Compute Jacobian matrix and determinant gradient with respect to each node coordinate
  *
  * dJ[l = i * 3 + j](3,3) = dN_i(3,1) * dxy_i(1,3)/dx_l (for node i and dimension j)
  */
-void Tetra4::buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eigen::MatrixXd> &dJ) const
+void Tetra4::buildGradientJ()
 {
-    MemTetra4 &mem = getMem();
-    Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
-
-    // Jacobian
-    dJ.resize(4 * 3); // 4 nodes and 3 dimensions (x, y, z)
-    for (size_t i = 0; i < 4; ++i)
+    CacheTetra4 &cache = getCache(order);
+    for (size_t k = 0; k < cache.getVGauss().getN(); ++k)
     {
-        for (size_t j = 0; j < 3; ++j)
+        Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
+
+        // Jacobian
+        gradJs[k].resize(4 * 3); // 4 nodes and 3 dimensions (x, y, z)
+        for (size_t i = 0; i < 4; ++i)
         {
-            dJ[i * 3 + j] = Eigen::Matrix3d::Zero();
-            dJ[i * 3 + j].col(j) = dsf.col(i); // retain sf gradients associated to ith node for dimension j
+            for (size_t j = 0; j < 3; ++j)
+            {
+                gradJs[k][i * 3 + j] = Eigen::Matrix3d::Zero();
+                gradJs[k][i * 3 + j].col(j) = dsf.col(i); // retain sf gradients associated to ith node for dimension j
+            }
         }
+        // Determinant
+        double detJ = this->getDetJ(k);
+        Eigen::MatrixXd const &iJ = this->getJinv(k);
+        gradDetJs[k].resize(4 * 3);
+        for (size_t i = 0; i < 4 * 3; ++i)
+            gradDetJs[k][i] = detJ * (iJ * gradJs[k][i]).trace();
     }
-    // Determinant
-    double detJ = mem.getDetJ(k);
-    Eigen::MatrixXd const &iJ = mem.getJinv(k);
-    dDetJ.resize(4 * 3);
-    for (size_t i = 0; i < dDetJ.size(); ++i)
-        dDetJ[i] = detJ * (iJ * dJ[i]).trace();
 }
 
 /**
  * @brief Return the element volume
  */
-double Tetra4::computeV() const
+void Tetra4::computeV()
 {
     CacheTetra4 &cache = getCache(order);
     GaussTetra4 &gauss = cache.gauss;
-    MemTetra4 &mem = getMem();
 
     // Gauss integration
-    double V = 0.0;
+    vol = 0.0;
     for (size_t k = 0; k < gauss.getN(); ++k)
-        V += gauss.getW(k) * mem.getDetJ(k);
-    return V;
+        vol += gauss.getW(k) * this->getDetJ(k);
 }
 /**
  * @brief Compute volume gradient with respect to each node coordinate
  *
  * dV[i] = w_k * dDetJ[i]_k
  */
-std::vector<double> Tetra4::buildGradientV(int dim) const
+void Tetra4::buildGradientV(int dim)
 {
     GaussTetra4 &gauss = getCache(order).gauss;
-    MemTetra4 &mem = getMem();
 
     // Gradient of volume
-    std::vector<double> dVol(4 * dim, 0.); // 4 nodes
+    gradVol.resize(4 * dim, 0.); // 4 nodes
     for (size_t k = 0; k < gauss.getN(); ++k)
     {
-        double w = gauss.getW(k);                       // GP weight
-        std::vector<double> dDetJ = mem.getGradDetJ(k); // gradient of Jacobian determinant
+        double w = gauss.getW(k);                         // GP weight
+        std::vector<double> dDetJ = this->getGradDetJ(k); // gradient of Jacobian determinant
         for (size_t i = 0; i < dDetJ.size(); ++i)
-            dVol[i] += w * dDetJ[i];
+            gradVol[i] += w * dDetJ[i];
     }
-    return dVol;
 }
 
 /**
@@ -143,5 +191,5 @@ Eigen::VectorXd Tetra4::computeGradient(std::vector<double> const &u, size_t k)
     for (size_t i = 0; i < nodes.size(); ++i)
         ue(i) = u[nodes[i]->row];
     // gradient in global axis: inv(J)_ij * d_jN_i * ue_i
-    return getMem().getJinv(k) * getCache(order).getDsf(k) * ue;
+    return this->getJinv(k) * getCache(order).getDsf(k) * ue;
 }
diff --git a/tbox/src/wTetra4.h b/tbox/src/wTetra4.h
index cc2bdb3970b4939322dfad8537dae06162d94798..371dcba8ae2b81b546be302805081cc297c253f0 100644
--- a/tbox/src/wTetra4.h
+++ b/tbox/src/wTetra4.h
@@ -32,30 +32,26 @@ namespace tbox
 class TBOX_API Tetra4 : public Element
 {
 #ifndef SWIG
-    friend class MemTetra4;
-
 private:
-    MemTetra4 *pmem; ///< private memory of precalculated values
     inline static CacheTetra4 &getCache(int n);
-    inline MemTetra4 &getMem() const;
 
 protected:
-    virtual double buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const override;
-    virtual void buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eigen::MatrixXd> &dJ) const override;
-
-public: // @todo
-    virtual double computeV() const override;
-    virtual std::vector<double> buildGradientV(int dim) const override;
+    virtual void buildJ() override;
+    virtual void buildGradientJ() override;
+    virtual void computeV() override;
+    virtual void buildGradientV(int dim) override;
 #endif
 
 public:
     Tetra4(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
-    virtual ~Tetra4();
+    virtual ~Tetra4() {}
     virtual ELTYPE type() const override { return ELTYPE::TETRA4; }
 
 #ifndef SWIG
+    virtual void initValues(bool isvol, int ordr = 2) override;
+    virtual void initGradients() override;
+    virtual void update() override;
     virtual Cache &getVCache() const override;
-    virtual Mem &getVMem() const override;
     virtual Eigen::VectorXd computeGradient(std::vector<double> const &u, size_t k) const override;
 #endif
 };
diff --git a/tbox/src/wTetra4.inl b/tbox/src/wTetra4.inl
index 506df073687d813473395d1e53cb172b695fb7b6..64d5d52a02cd64bf155669a5f609ee955e1851dc 100644
--- a/tbox/src/wTetra4.inl
+++ b/tbox/src/wTetra4.inl
@@ -14,14 +14,6 @@
  * limitations under the License.
  */
 
-/**
- * @brief Return private memory for Tetra4 private methods
- */
-inline MemTetra4 &Tetra4::getMem() const
-{
-    return *pmem;
-}
-
 /**
  * @brief Return cache for Tetra4 private methods
  */
diff --git a/tbox/src/wTri3.cpp b/tbox/src/wTri3.cpp
index f145c4ed2bc8b947ed188c6e1c990f7f54e34849..ea16dd09ad66fc03aff6711d6cc5cd8fdd12db5a 100644
--- a/tbox/src/wTri3.cpp
+++ b/tbox/src/wTri3.cpp
@@ -18,7 +18,6 @@
 #include "wSfTri3.h"
 #include "wGaussTri3.h"
 #include "wCacheTri3.h"
-#include "wMemTri3.h"
 #include "wNode.h"
 
 using namespace tbox;
@@ -26,147 +25,214 @@ using namespace tbox;
 Tri3::Tri3(int n, Tag *_ptag, Tag *_etag,
            std::vector<Tag *> &_parts, std::vector<Node *> &nods) : Element(n, _ptag, _etag, _parts, nods)
 {
-    pmem = new MemTri3(*this, static_cast<int>(getCache(order).gauss.getN()));
 }
-Tri3::~Tri3()
+
+/**
+ * @brief Return cache holding shape functions and Gauss points for all Tri3 elements
+ */
+Cache &Tri3::getVCache() const
 {
-    if (pmem)
-        delete pmem;
+    return getCache(order);
 }
 
 /**
- * @brief Return private memory handling Jacobian for this Tri3 element
+ * @brief Initialize precomputed values
  */
-Mem &Tri3::getVMem() const
+void Tri3::initValues(bool isvol, int ordr)
 {
-    return getMem();
+    // Set flags
+    isVol = isvol;
+    order = ordr;
+    hasVals = true;
+    // Update values
+    size_t ngp = getCache(order).getVGauss().getN();
+    Js.resize(ngp);
+    iJs.resize(ngp);
+    detJs.resize(ngp);
+    this->update();
 }
-
 /**
- * @brief Return cache holding shape functions and Gauss points for all Tri3 elements
+ * @brief Initialize precomputed gradients
  */
-Cache &Tri3::getVCache() const
+void Tri3::initGradients()
 {
-    return getCache(order);
+    // Set flags
+    hasGrads = true;
+    // Update values
+    size_t ngp = getCache(order).getVGauss().getN();
+    gradJs.resize(ngp);
+    gradDetJs.resize(ngp);
+    this->update();
+}
+/**
+ * @brief Update precomputed values and gradients
+ */
+void Tri3::update()
+{
+    // Update values
+    if (hasVals)
+    {
+        // Jacobian and normal
+        // Tri3 is not of same dimension as Problem (surface element)
+        if (!isVol)
+        {
+            this->buildSurfaceJ();
+            this->computeNormal();
+        }
+        // Tri3 is of same dimension as Problem (volume element)
+        else
+            this->buildJ();
+        // Volume and CG
+        this->computeV();
+        this->computeCg();
+    }
+    // Update gradients
+    if (hasGrads)
+    {
+        // Jacobian, volume and normal
+        // Tri3 is not of same dimension as Problem (surface element)
+        if (!isVol)
+        {
+            this->buildGradientSurfaceJ();
+            this->buildGradientV(3);
+            this->computeGradientNormal();
+        }
+        // Tri3 is of same dimension as Problem (volume element)
+        else
+        {
+            this->buildGradientJ();
+            this->buildGradientV(2);
+        }
+    }
 }
 
 /**
- * @brief Compute Jacobian determinant at Gauss point k
+ * @brief Compute Jacobian determinant
  *
  * detJ = norm(x_i dN_i,xi X x_i dN_i,eta)
  */
-double Tri3::buildJ(size_t k) const
+void Tri3::buildSurfaceJ()
 {
-    std::vector<Eigen::Vector3d> ts = this->buildTangents(k);
-    return ts[0].cross(ts[1]).norm();
+    for (size_t k = 0; k < getCache(order).getVGauss().getN(); ++k)
+    {
+        std::vector<Eigen::Vector3d> ts = this->buildTangents(k);
+        detJs[k] = ts[0].cross(ts[1]).norm();
+    }
 }
 /**
- * @brief Compute Jacobian matrix, inverse and determinant at Gauss point k
+ * @brief Compute Jacobian matrix, inverse and determinant
  *
  * J_ij(2,2) = diN_n * xj_n (n is the node)
  */
-double Tri3::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
+void Tri3::buildJ()
 {
-    Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
+    CacheTri3 &cache = getCache(order);
+    for (size_t k = 0; k < cache.getVGauss().getN(); ++k)
+    {
+        Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
 
-    Eigen::Matrix2d JJ = Eigen::Matrix2d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
-    size_t i = 0;
-    for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
-        JJ += dff.col(i) * (*it)->pos.transpose().block<1, 2>(0, 0);
-    J = JJ;
-    iJ = JJ.inverse();
-    return JJ.determinant();
+        Eigen::Matrix2d JJ = Eigen::Matrix2d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
+        size_t i = 0;
+        for (auto it = nodes.begin(); it != nodes.end(); ++it, ++i)
+            JJ += dff.col(i) * (*it)->pos.transpose().block<1, 2>(0, 0);
+        Js[k] = JJ;
+        iJs[k] = JJ.inverse();
+        detJs[k] = JJ.determinant();
+        this->checkJac(k);
+    }
 }
 
 /**
- * @brief Compute Jacobian determinant gradient with respect to each node coordinate, at Gauss point k
+ * @brief Compute Jacobian determinant gradient with respect to each node coordinate
  *
  * dDetJ[l = i * 3 + j] = detJ_l * (dt0 \cross t1 + t0 \cross dt1)_l;
  */
-void Tri3::buildGradientJ(size_t k, std::vector<double> &dDetJ) const
+void Tri3::buildGradientSurfaceJ()
 {
-    Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
-
-    // Tangent vectors and determinant
-    std::vector<Eigen::Vector3d> ts = this->buildTangents(k);
-    Eigen::Vector3d detJ = ts[0].cross(ts[1]).normalized();
-
-    // Jacobian gradient
-    dDetJ.resize(3 * 3); // 3 nodes and 3 dimensions (x, y, z)
-    for (size_t i = 0; i < 3; ++i)
+    CacheTri3 &cache = getCache(order);
+    for (size_t k = 0; k < getCache(order).getVGauss().getN(); ++k)
     {
-        for (size_t j = 0; j < 3; ++j)
+        // Shape functions, tangent vectors and determinant
+        Eigen::MatrixXd const &dsf = cache.getDsf(k);
+        std::vector<Eigen::Vector3d> ts = this->buildTangents(k);
+        Eigen::Vector3d detJ = ts[0].cross(ts[1]).normalized();
+
+        // Jacobian gradient
+        gradDetJs[k].resize(3 * 3); // 3 nodes and 3 dimensions (x, y, z)
+        for (size_t i = 0; i < 3; ++i)
         {
-            Eigen::Vector3d dt0 = Eigen::Vector3d::Zero(), dt1 = Eigen::Vector3d::Zero();
-            dt0(j) = dsf(0, i);
-            dt1(j) = dsf(1, i); // retain sf gradients associated to ith node for dimension j
-            dDetJ[i * 3 + j] = detJ.dot(dt0.cross(ts[1]) + ts[0].cross(dt1));
+            for (size_t j = 0; j < 3; ++j)
+            {
+                Eigen::Vector3d dt0 = Eigen::Vector3d::Zero(), dt1 = Eigen::Vector3d::Zero();
+                dt0(j) = dsf(0, i);
+                dt1(j) = dsf(1, i); // retain sf gradients associated to ith node for dimension j
+                gradDetJs[k][i * 3 + j] = detJ.dot(dt0.cross(ts[1]) + ts[0].cross(dt1));
+            }
         }
     }
 }
 /**
- * @brief Compute Jacobian matrix gradient with respect to each node coordinate, at Gauss point k
+ * @brief Compute Jacobian matrix gradient with respect to each node coordinate
  *
  * dJ[l = i * 2 + j](2,2) = dN_i(2,1) * dxy_i(1,2)/dx_l (for node i and dimension j)
  */
-void Tri3::buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eigen::MatrixXd> &dJ) const
+void Tri3::buildGradientJ()
 {
-    MemTri3 &mem = getMem();
-    Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
-
-    // Jacobian
-    dJ.resize(3 * 2); // 3 nodes and 2 dimensions (x, y)
-    for (size_t i = 0; i < 3; ++i)
+    CacheTri3 &cache = getCache(order);
+    for (size_t k = 0; k < cache.getVGauss().getN(); ++k)
     {
-        for (size_t j = 0; j < 2; ++j)
+        Eigen::MatrixXd const &dsf = cache.getDsf(k);
+
+        // Jacobian
+        gradJs[k].resize(3 * 2); // 3 nodes and 2 dimensions (x, y)
+        for (size_t i = 0; i < 3; ++i)
         {
-            dJ[i * 2 + j] = Eigen::Matrix2d::Zero();
-            dJ[i * 2 + j].col(j) = dsf.col(i); // retain sf gradients associated to ith node for dimension j
+            for (size_t j = 0; j < 2; ++j)
+            {
+                gradJs[k][i * 2 + j] = Eigen::Matrix2d::Zero();
+                gradJs[k][i * 2 + j].col(j) = dsf.col(i); // retain sf gradients associated to ith node for dimension j
+            }
         }
+        // Determinant
+        double detJ = this->getDetJ(k);
+        Eigen::MatrixXd const &iJ = this->getJinv(k);
+        gradDetJs[k].resize(3 * 2);
+        for (size_t i = 0; i < 3 * 2; ++i)
+            gradDetJs[k][i] = detJ * (iJ * gradJs[k][i]).trace();
     }
-    // Determinant
-    double detJ = mem.getDetJ(k);
-    Eigen::MatrixXd const &iJ = mem.getJinv(k);
-    dDetJ.resize(3 * 2);
-    for (size_t i = 0; i < dDetJ.size(); ++i)
-        dDetJ[i] = detJ * (iJ * dJ[i]).trace();
 }
 
 /**
  * @brief Return the element surface
  */
-double Tri3::computeV() const
+void Tri3::computeV()
 {
     CacheTri3 &cache = getCache(order);
     GaussTri3 &gauss = cache.gauss;
-    MemTri3 &mem = getMem();
 
     // Gauss integration
-    double V = 0.0;
+    vol = 0.0;
     for (size_t k = 0; k < gauss.getN(); ++k)
-        V += gauss.getW(k) * mem.getDetJ(k);
-    return V;
+        vol += gauss.getW(k) * this->getDetJ(k);
 }
 /**
  * @brief Compute volume gradient with respect to each node coordinate
  *
  * dV[i] = w_k * dDetJ[i]_k
  */
-std::vector<double> Tri3::buildGradientV(int dim) const
+void Tri3::buildGradientV(int dim)
 {
     GaussTri3 &gauss = getCache(order).gauss;
-    MemTri3 &mem = getMem();
 
     // Gradient of volume
-    std::vector<double> dVol(3 * dim, 0.); // 3 nodes
+    gradVol.resize(3 * dim, 0.); // 3 nodes
     for (size_t k = 0; k < gauss.getN(); ++k)
     {
-        double w = gauss.getW(k);                       // GP weight
-        std::vector<double> dDetJ = mem.getGradDetJ(k); // gradient of Jacobian determinant
+        double w = gauss.getW(k);                         // GP weight
+        std::vector<double> dDetJ = this->getGradDetJ(k); // gradient of Jacobian determinant
         for (size_t i = 0; i < dDetJ.size(); ++i)
-            dVol[i] += w * dDetJ[i];
+            gradVol[i] += w * dDetJ[i];
     }
-    return dVol;
 }
 
 /**
@@ -201,16 +267,16 @@ std::vector<Eigen::Vector3d> Tri3::buildTangents(size_t k) const
 /**
  * @brief Compute element unit normal vector
  */
-Eigen::Vector3d Tri3::normal() const
+void Tri3::computeNormal()
 {
     std::vector<Eigen::Vector3d> ts = this->computeTangents();
-    return ts[0].cross(ts[1]).normalized();
+    normal = ts[0].cross(ts[1]).normalized();
 }
 
 /**
  * @brief Compute gradient of element unit normal vector with respect to each node coordinate
  */
-std::vector<Eigen::Vector3d> Tri3::gradientNormal() const
+void Tri3::computeGradientNormal()
 {
     // Compute normal and norm
     std::vector<Eigen::Vector3d> ts = this->computeTangents();
@@ -219,17 +285,16 @@ std::vector<Eigen::Vector3d> Tri3::gradientNormal() const
 
     // Compute gradients
     Eigen::Matrix3d inn = (Eigen::Matrix3d::Identity() - n * n.transpose() / (nn * nn)) / nn; // gradient of n/|n|
-    std::vector<Eigen::Vector3d> dn(3 * 3);                                                   // 3 nodes and 3 dimensions
+    gradNormal.resize(3 * 3);                                                                 // 3 nodes and 3 dimensions
     for (size_t i = 0; i < 3; ++i)
     {
         for (size_t j = 0; j < 3; ++j)
         {
             Eigen::Matrix3d dx = Eigen::Matrix3d::Zero();
-            dx(j, i) = 1;                                                                                        // dpos_i / dx_j
-            dn[i * 3 + j] = inn * ((dx.col(1) - dx.col(0)).cross(ts[1]) + (ts[0]).cross(dx.col(2) - dx.col(0))); // outer gradient times gradient of cross product
+            dx(j, i) = 1;                                                                                                // dpos_i / dx_j
+            gradNormal[i * 3 + j] = inn * ((dx.col(1) - dx.col(0)).cross(ts[1]) + (ts[0]).cross(dx.col(2) - dx.col(0))); // outer gradient times gradient of cross product
         }
     }
-    return dn;
 }
 
 /**
@@ -242,5 +307,5 @@ Eigen::VectorXd Tri3::computeGradient(std::vector<double> const &u, size_t k) co
     for (size_t i = 0; i < nodes.size(); ++i)
         ue(i) = u[nodes[i]->row];
     // Gradient in global axis: inv(J)_ij * d_jN_i * ue_i
-    return getMem().getJinv(k) * getCache(order).getDsf(k) * ue;
+    return this->getJinv(k) * getCache(order).getDsf(k) * ue;
 }
diff --git a/tbox/src/wTri3.h b/tbox/src/wTri3.h
index bb82eab1bb13448d34c42172d74eb1476c92a0d9..c042c4c5a93094652fadab0e549d35c29fb4c9c3 100644
--- a/tbox/src/wTri3.h
+++ b/tbox/src/wTri3.h
@@ -32,35 +32,31 @@ namespace tbox
 class TBOX_API Tri3 : public Element
 {
 #ifndef SWIG
-    friend class MemTri3;
-
 private:
-    MemTri3 *pmem; ///< private memory of precalculated values
-    inline MemTri3 &getMem() const;
     inline static CacheTri3 &getCache(int order);
 
 protected:
-    virtual double buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const override;
-    virtual double buildJ(size_t k) const override;
-    virtual void buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eigen::MatrixXd> &dJ) const override;
-    virtual void buildGradientJ(size_t k, std::vector<double> &dDetJ) const override;
-
-public: // @todo
-    virtual double computeV() const override;
-    virtual std::vector<double> buildGradientV(int dim) const override;
-    virtual Eigen::Vector3d normal() const override;
-    virtual std::vector<Eigen::Vector3d> gradientNormal() const override;
+    virtual void buildJ() override;
+    virtual void buildSurfaceJ() override;
+    virtual void buildGradientJ() override;
+    virtual void buildGradientSurfaceJ() override;
+    virtual void computeV() override;
+    virtual void buildGradientV(int dim) override;
+    virtual void computeNormal() override;
+    virtual void computeGradientNormal() override;
     virtual std::vector<Eigen::Vector3d> computeTangents() const override;
     virtual std::vector<Eigen::Vector3d> buildTangents(size_t k) const override;
 #endif
 
 public:
     Tri3(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
-    virtual ~Tri3();
+    virtual ~Tri3() {}
     virtual ELTYPE type() const override { return ELTYPE::TRI3; }
 
 #ifndef SWIG
-    virtual Mem &getVMem() const override;
+    virtual void initValues(bool isvol, int ordr = 2) override;
+    virtual void initGradients() override;
+    virtual void update() override;
     virtual Cache &getVCache() const override;
     virtual Eigen::VectorXd computeGradient(std::vector<double> const &u, size_t k) const override;
 #endif
diff --git a/tbox/src/wTri3.inl b/tbox/src/wTri3.inl
index 23df4f1e2030610350cc0d48008591d08540dd2f..57f12736e0ee17678604e876113fd524e4e4de8c 100644
--- a/tbox/src/wTri3.inl
+++ b/tbox/src/wTri3.inl
@@ -14,14 +14,6 @@
  * limitations under the License.
  */
 
-/**
- * @brief Return private memory for Tri3 private methods
- */
-inline MemTri3 &Tri3::getMem() const
-{
-    return *pmem;
-}
-
 /**
  * @brief Return cache for Tri3 private methods
  */
diff --git a/tbox/src/wUnitTest.cpp b/tbox/src/wUnitTest.cpp
index e3f50dae83a085af3e91ab531d75f9ffc9cedc14..7b886c29e77e50893e9d5c64991e1b34e538b72d 100644
--- a/tbox/src/wUnitTest.cpp
+++ b/tbox/src/wUnitTest.cpp
@@ -22,7 +22,6 @@
 #include "wGmshExport.h"
 #include "wResults.h"
 #include "wMshDeform.h"
-#include "wMem.h"
 #include <Eigen/Dense>
 
 using namespace tbox;
@@ -84,7 +83,7 @@ void UnitTest::mshDeform(MshDeform *mshDef, int d)
     // Initialize memory (volume element only)
     for (auto e : mshDef->msh->elems)
         if ((d == 2 && e->type() == ELTYPE::TRI3) || (d == 3 && e->type() == ELTYPE::TETRA4))
-            e->getVMem().update(true);
+            e->initValues(true);
     // Call deform
     mshDef->deform();
 }
diff --git a/waves/src/wProblem.cpp b/waves/src/wProblem.cpp
index b7a4a891e742a457cd817655bffbbfcb38907701..82907475095639af3f19761807974136ba718222 100644
--- a/waves/src/wProblem.cpp
+++ b/waves/src/wProblem.cpp
@@ -20,7 +20,7 @@
 #include "wMedium.h"
 #include "wBoundary.h"
 #include "wElement.h"
-#include "wMem.h"
+
 #include "wTag.h"
 using namespace waves;
 
@@ -37,19 +37,18 @@ Problem::~Problem()
 }
 
 /**
- * @brief Update the elements memory (Jacobian)
- * @authors Adrien Crovato
+ * @brief Initialize the elements precomputed values
  */
-void Problem::updateMem()
+void Problem::initElems()
 {
     // Update volume Jacobian
     for (auto vol : media)
         for (auto e : vol->tag->elems)
-            e->getVMem().update(true);
+            e->initValues(true);
     // Update surface Jacobian (Neumann B.C.)
     for (auto surf : bnds)
         for (auto e : surf->tag->elems)
-            e->getVMem().update(false);
+            e->initValues(false);
 }
 
 void Problem::write(std::ostream &out) const
diff --git a/waves/src/wProblem.h b/waves/src/wProblem.h
index aaf5b9637132b0f64c9534903d6eadb9b1262858..0f9e9b04f7543086a614ab087fdf0ace9d4642d6 100644
--- a/waves/src/wProblem.h
+++ b/waves/src/wProblem.h
@@ -49,7 +49,7 @@ public:
     void add(std::shared_ptr<Boundary> b) { bnds.push_back(b); }
 
 #ifndef SWIG
-    void updateMem();
+    void initElems();
     virtual void write(std::ostream &out) const override;
 #endif
 };
diff --git a/waves/src/wTimeIntegration.cpp b/waves/src/wTimeIntegration.cpp
index d958627397747f8a6bd35e424fcfcee28f11d3ed..9d3c476cd32df70bd60cb1ffdbd650d0c3faaa8b 100644
--- a/waves/src/wTimeIntegration.cpp
+++ b/waves/src/wTimeIntegration.cpp
@@ -52,7 +52,7 @@ TimeIntegration::TimeIntegration(std::shared_ptr<Problem> _pbl) : pbl(_pbl)
     stopit = false;
 
     // Update element memory
-    pbl->updateMem();
+    pbl->initElems();
 }
 
 void TimeIntegration::start(MshExport *mshWriter)
@@ -101,7 +101,7 @@ void TimeIntegration::buildS(Eigen::SparseMatrix<double, Eigen::RowMajor> &S2)
     for (auto bnd : pbl->bnds)
     {
         std::cout << "\tprocessing " << *bnd << '\n';
-        std::cout << "normal=" << static_cast<Quad4 *>(bnd->tag->elems[0])->normal() << '\n';
+        std::cout << "normal=" << static_cast<Quad4 *>(bnd->tag->elems[0])->getNormal() << '\n';
         tbb::parallel_do(bnd->tag->elems.begin(), bnd->tag->elems.end(),
                          [&](Element *e) {
                              if (e->type() != ELTYPE::QUAD4)
diff --git a/waves/src/wWaveTerm.cpp b/waves/src/wWaveTerm.cpp
index e39582ec569b9782b05a1324817b0b2f6de3306e..8cb036fc86a8259809491a7182a3ca1473fc58a6 100644
--- a/waves/src/wWaveTerm.cpp
+++ b/waves/src/wWaveTerm.cpp
@@ -19,7 +19,6 @@
 #include "wElement.h"
 #include "wCache.h"
 #include "wGauss.h"
-#include "wMem.h"
 
 using namespace waves;
 using namespace tbox;
@@ -32,18 +31,17 @@ Eigen::MatrixXd WaveTerm::buildK(Element const &e, std::vector<double> const &u)
     // Get precomputed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Stiffness matrix
     Eigen::MatrixXd K = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
     for (size_t k = 0; k < gauss.getN(); ++k)
     {
         // Jacobian inverse and shape functions
-        Eigen::Matrix3d const &J = mem.getJinv(k);
+        Eigen::Matrix3d const &J = e.getJinv(k);
         Eigen::MatrixXd const &dff = cache.getDsf(k);
 
         // Elementary stiffness matrix
-        K += (J * dff).transpose() * J * dff * gauss.getW(k) * mem.getDetJ(k);
+        K += (J * dff).transpose() * J * dff * gauss.getW(k) * e.getDetJ(k);
     }
     return K;
 }
@@ -56,7 +54,6 @@ Eigen::MatrixXd WaveTerm::buildM(Element const &e)
     // Get precomputed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    Mem &mem = e.getVMem();
 
     // Mass matrix
     Eigen::MatrixXd M = Eigen::MatrixXd::Zero(e.nodes.size(), e.nodes.size());
@@ -64,7 +61,7 @@ Eigen::MatrixXd WaveTerm::buildM(Element const &e)
     {
         // Shape functions
         Eigen::VectorXd const &ff = cache.getSf(k);
-        M += ff * ff.transpose() * gauss.getW(k) * mem.getDetJ(k);
+        M += ff * ff.transpose() * gauss.getW(k) * e.getDetJ(k);
     }
     return M;
 }