diff --git a/tbox/src/wCache.h b/tbox/src/wCache.h
index e5dab88ff7edd9d8d7a644093ea96a95c5b9a20e..f20921ce6c6e9d7013dbefea2bddc297c8c434f7 100644
--- a/tbox/src/wCache.h
+++ b/tbox/src/wCache.h
@@ -37,7 +37,7 @@ protected:
     std::vector<Eigen::VectorXd> ff;  ///< shape functions
 
 public:
-    explicit Cache(int _npg);
+    explicit Cache(int _n);
     virtual ~Cache() {}
 
     inline Eigen::MatrixXd const &getDsf(size_t k);
diff --git a/tbox/src/wCacheHex8.h b/tbox/src/wCacheHex8.h
index b0f6e9a779190965967465e90453b3e01a7756ff..212dba307b2dbc29cf9855dd871a01f35c9af6b8 100644
--- a/tbox/src/wCacheHex8.h
+++ b/tbox/src/wCacheHex8.h
@@ -36,7 +36,7 @@ public:
     SfHex8 sf;
     GaussHex8 gauss;
 
-    explicit CacheHex8(int _npg = 2);
+    explicit CacheHex8(int _n = 2);
 
     virtual Gauss &getVGauss() override;
 };
diff --git a/tbox/src/wCacheLine2.h b/tbox/src/wCacheLine2.h
index 1d0f160fb2919d860b27b0d8cb838be50307a767..a912013a09f760f6114835a05e2241c1de55e3b7 100644
--- a/tbox/src/wCacheLine2.h
+++ b/tbox/src/wCacheLine2.h
@@ -36,7 +36,7 @@ public:
     SfLine2 sf;
     GaussLine2 gauss;
 
-    explicit CacheLine2(int _npg = 2);
+    explicit CacheLine2(int _n = 2);
 
     virtual Gauss &getVGauss() override;
 };
diff --git a/tbox/src/wCachePoint1.h b/tbox/src/wCachePoint1.h
index d14baf726376211e521ef5233a21c662278e8468..9d5f290b1f2c46d89e8d3c926d52eab2eb0d45d3 100644
--- a/tbox/src/wCachePoint1.h
+++ b/tbox/src/wCachePoint1.h
@@ -33,7 +33,7 @@ namespace tbox
 class TBOX_API CachePoint1 : public Cache
 {
 public:
-    explicit CachePoint1(int _npg = 0);
+    explicit CachePoint1(int _n = 0);
 };
 
 } // namespace tbox
diff --git a/tbox/src/wCacheQuad4.h b/tbox/src/wCacheQuad4.h
index 11cc0210502120d5a68d6c4974eeacd863912017..49a2fb1e5ca9d38e6ac851f02aa957389240a7d1 100644
--- a/tbox/src/wCacheQuad4.h
+++ b/tbox/src/wCacheQuad4.h
@@ -36,7 +36,7 @@ public:
     SfQuad4 sf;
     GaussQuad4 gauss;
 
-    explicit CacheQuad4(int _npg = 2);
+    explicit CacheQuad4(int _n = 2);
 
     virtual Gauss &getVGauss() override;
 };
diff --git a/tbox/src/wCacheTetra4.h b/tbox/src/wCacheTetra4.h
index b56b1b1b95a503311401632f36089e8e428cbe14..4d1062ff11d2f0d92b0022afbad1da79247b4151 100644
--- a/tbox/src/wCacheTetra4.h
+++ b/tbox/src/wCacheTetra4.h
@@ -36,7 +36,7 @@ public:
     SfTetra4 sf;
     GaussTetra4 gauss;
 
-    explicit CacheTetra4(int _npg = 4);
+    explicit CacheTetra4(int _n = 2);
 
     virtual Gauss &getVGauss() override;
 };
diff --git a/tbox/src/wCacheTri3.h b/tbox/src/wCacheTri3.h
index da73a5d72f3a7a34a315fdaa4a055b977c8ccd56..c7491eb7bf8822d4555337c47e6b0e11a22761f5 100644
--- a/tbox/src/wCacheTri3.h
+++ b/tbox/src/wCacheTri3.h
@@ -36,7 +36,7 @@ public:
     SfTri3 sf;
     GaussTri3 gauss;
 
-    explicit CacheTri3(int _npg = 3);
+    explicit CacheTri3(int _n = 2);
 
     virtual Gauss &getVGauss() override;
 };
diff --git a/tbox/src/wElement.cpp b/tbox/src/wElement.cpp
index c3cd11273dec76ba33058613a17de2aedf769140..a1b3d5615ba43ad9df9d2dec4ed60ac5a0e8adfb 100644
--- a/tbox/src/wElement.cpp
+++ b/tbox/src/wElement.cpp
@@ -84,7 +84,8 @@ operator<<(std::ostream &out, MATTYPE const &obj)
 
 Element::Element(int n, Tag *_ptag, Tag *_etag,
                  std::vector<Tag *> &_parts, std::vector<Node *> &nods)
-    : wObject(), no(n), ptag(_ptag), etag(_etag), parts(_parts), nodes(nods)
+    : wObject(), order(2), no(n), ptag(_ptag), etag(_etag), parts(_parts),
+      nodes(nods)
 {
     if (ptag)
         ptag->elems.push_back(this);
diff --git a/tbox/src/wElement.h b/tbox/src/wElement.h
index 250884efc9da5f8cfebf0b9c424f70f858a8468c..48d60c9252ac06f959b6ff80de8d3220c3b876e5 100644
--- a/tbox/src/wElement.h
+++ b/tbox/src/wElement.h
@@ -59,6 +59,8 @@ TBOX_API std::ostream &operator<<(std::ostream &out, MATTYPE const &obj);
  */
 class TBOX_API Element : public fwk::wObject
 {
+protected:
+    int order; ///< order of integration (Gauss quadrature)
 public:
     int no;                    ///< label
     Tag *ptag;                 ///< physical tag   (group)
diff --git a/tbox/src/wGaussTetra4.cpp b/tbox/src/wGaussTetra4.cpp
index 5bdff84d13697510982c64966c6442523b9f08a2..0e1a4e351492319a0cbe20869132748c6adcf845 100644
--- a/tbox/src/wGaussTetra4.cpp
+++ b/tbox/src/wGaussTetra4.cpp
@@ -24,6 +24,15 @@ GaussTetra4::GaussTetra4(int n) : Gauss()
 {
     switch (n)
     {
+    case 1:
+    {
+        p.resize(1);
+        w.resize(1);
+        ngp = 1;
+        p[0] = Eigen::Vector3d(0.25, 0.25, 0.25);
+        w[0] = 1. / 6.;
+        break;
+    }
     case 2:
     {
         p.resize(4);
diff --git a/tbox/src/wHex8.cpp b/tbox/src/wHex8.cpp
index 792c5647b28d0345e89deffb134ba96f545a205d..da1332f572e05d3ca27238dc4604c04b581ce29a 100644
--- a/tbox/src/wHex8.cpp
+++ b/tbox/src/wHex8.cpp
@@ -27,7 +27,7 @@ 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().gauss.getN()));
+    pmem = new MemHex8(*this, static_cast<int>(getCache(order).gauss.getN()));
 }
 Hex8::~Hex8()
 {
@@ -48,7 +48,7 @@ Mem &Hex8::getVMem() const
  */
 Cache &Hex8::getVCache() const
 {
-    return getCache();
+    return getCache(order);
 }
 
 /**
@@ -58,7 +58,7 @@ Cache &Hex8::getVCache() const
  */
 double Hex8::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
 {
-    Eigen::MatrixXd const &dff = getCache().getDsf(k);
+    Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
 
     Eigen::Matrix3d JJ = Eigen::Matrix3d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
     size_t i = 0;
@@ -76,7 +76,7 @@ double Hex8::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
  */
 Eigen::MatrixXd Hex8::buildK(std::vector<double> const &u, Fct0 const &fct)
 {
-    CacheHex8 &cache = getCache();
+    CacheHex8 &cache = getCache(order);
     GaussHex8 &gauss = cache.gauss;
     MemHex8 &mem = getMem();
 
@@ -101,7 +101,7 @@ Eigen::MatrixXd Hex8::buildK(std::vector<double> const &u, Fct0 const &fct)
  */
 Eigen::MatrixXd Hex8::buildK_mechanical(Eigen::MatrixXd const &H)
 {
-    CacheHex8 &cache = getCache();
+    CacheHex8 &cache = getCache(order);
     GaussHex8 &gauss = cache.gauss;
     MemHex8 &mem = getMem();
 
@@ -139,7 +139,7 @@ Eigen::MatrixXd Hex8::buildK_mechanical(Eigen::MatrixXd const &H)
  */
 Eigen::MatrixXd Hex8::buildS_thermomechanical(Eigen::MatrixXd const &D)
 {
-    CacheHex8 &cache = getCache();
+    CacheHex8 &cache = getCache(order);
     GaussHex8 &gauss = cache.gauss;
     MemHex8 &mem = getMem();
 
@@ -163,7 +163,7 @@ Eigen::MatrixXd Hex8::buildS_thermomechanical(Eigen::MatrixXd const &D)
  */
 Eigen::MatrixXd Hex8::buildL_thermal(Eigen::MatrixXd const &C)
 {
-    CacheHex8 &cache = getCache();
+    CacheHex8 &cache = getCache(order);
     GaussHex8 &gauss = cache.gauss;
     MemHex8 &mem = getMem();
 
@@ -185,7 +185,7 @@ Eigen::MatrixXd Hex8::buildL_thermal(Eigen::MatrixXd const &C)
 
 void Hex8::strain_stress(Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, Eigen::MatrixXd const &H, std::vector<double> const &u)
 {
-    CacheHex8 &cache = getCache();
+    CacheHex8 &cache = getCache(order);
     GaussHex8 &gauss = cache.gauss;
     MemHex8 &mem = getMem();
 
@@ -222,7 +222,7 @@ void Hex8::strain_stress(Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, Eigen
  */
 Eigen::MatrixXd Hex8::buildM()
 {
-    CacheHex8 &cache = getCache();
+    CacheHex8 &cache = getCache(order);
     GaussHex8 &gauss = cache.gauss;
     MemHex8 &mem = getMem();
 
@@ -244,7 +244,7 @@ Eigen::MatrixXd Hex8::buildM()
  */
 double Hex8::computeV(std::vector<double> const &u, Fct0 const &fct)
 {
-    CacheHex8 &cache = getCache();
+    CacheHex8 &cache = getCache(order);
     GaussHex8 &gauss = cache.gauss;
     MemHex8 &mem = getMem();
 
diff --git a/tbox/src/wHex8.h b/tbox/src/wHex8.h
index 7850cc980c73aa08807fed8c14f97d2e6b0cfafe..1c563df033a1ca4381e156159434978d06f87461 100644
--- a/tbox/src/wHex8.h
+++ b/tbox/src/wHex8.h
@@ -34,7 +34,7 @@ class TBOX_API Hex8 : public Element
 #ifndef SWIG
     friend class MemHex8;
     MemHex8 *pmem; ///< private memory of precalculated values
-    inline static CacheHex8 &getCache();
+    inline static CacheHex8 &getCache(int n);
     inline MemHex8 &getMem() const;
     virtual double buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const override;
 
diff --git a/tbox/src/wHex8.inl b/tbox/src/wHex8.inl
index aa6e937f0d0d57093d9cfa69f92e23925bbe2780..6912e483f986152c6952c5bc7ef2c45372ebe2e5 100644
--- a/tbox/src/wHex8.inl
+++ b/tbox/src/wHex8.inl
@@ -25,8 +25,8 @@ inline MemHex8 &Hex8::getMem() const
 /**
  * @brief Return cache for Hex8 private methods
  */
-inline CacheHex8 &Hex8::getCache()
+inline CacheHex8 &Hex8::getCache(int n)
 {
-    static Lazy<CacheHex8> cache(2);
+    static Lazy<CacheHex8> cache(n);
     return cache.get();
 }
diff --git a/tbox/src/wLine2.cpp b/tbox/src/wLine2.cpp
index 00dff4adc0f39e4b4442e0674f891901b0c5c1a7..986904ef6cb4d0cf3584eac14d14173d87bc4ebe 100644
--- a/tbox/src/wLine2.cpp
+++ b/tbox/src/wLine2.cpp
@@ -28,7 +28,7 @@ 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().gauss.getN()));
+    pmem = new MemLine2(*this, static_cast<int>(getCache(order).gauss.getN()));
 }
 Line2::~Line2()
 {
@@ -49,7 +49,7 @@ Mem &Line2::getVMem() const
  */
 Cache &Line2::getVCache() const
 {
-    return getCache();
+    return getCache(order);
 }
 
 /**
@@ -65,7 +65,7 @@ std::vector<Eigen::Vector3d> Line2::computeTangents() const
  */
 std::vector<Eigen::Vector3d> Line2::buildTangents(size_t k) const
 {
-    Eigen::MatrixXd const &dff = getCache().getDsf(k);
+    Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
 
     std::vector<Eigen::Vector3d> t(1, Eigen::Vector3d::Zero());
     size_t i = 0;
@@ -122,7 +122,7 @@ double Line2::buildJ(size_t k) const
  */
 void Line2::buildGradientJ(size_t k, std::vector<double> &dDetJ) const
 {
-    Eigen::MatrixXd const &dsf = getCache().getDsf(k);
+    Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
 
     // Tangent vectors and determinant
     Eigen::Vector3d detJ = this->buildTangents(k)[0].normalized();
@@ -147,7 +147,7 @@ void Line2::buildGradientJ(size_t k, std::vector<double> &dDetJ) const
  */
 std::vector<double> Line2::buildGradientV(int dim) const
 {
-    GaussLine2 &gauss = getCache().gauss;
+    GaussLine2 &gauss = getCache(order).gauss;
     MemLine2 &mem = getMem();
 
     // Gradient of volume
@@ -169,7 +169,7 @@ std::vector<double> Line2::buildGradientV(int dim) const
  */
 Eigen::VectorXd Line2::builds(std::vector<double> const &u, Fct0 const &f)
 {
-    CacheLine2 &cache = getCache();
+    CacheLine2 &cache = getCache(order);
     GaussLine2 &gauss = cache.gauss;
     MemLine2 &mem = getMem();
 
@@ -187,7 +187,7 @@ Eigen::VectorXd Line2::builds(std::vector<double> const &u, Fct0 const &f)
  */
 Eigen::VectorXd Line2::builds(std::vector<double> const &u, Fct1 const &f)
 {
-    CacheLine2 &cache = getCache();
+    CacheLine2 &cache = getCache(order);
     GaussLine2 &gauss = cache.gauss;
     MemLine2 &mem = getMem();
 
@@ -216,7 +216,7 @@ double Line2::computeGrad2(std::vector<double> const &u, size_t k)
  */
 double Line2::computeV(std::vector<double> const &u, Fct0 const &fct)
 {
-    CacheLine2 &cache = getCache();
+    CacheLine2 &cache = getCache(order);
     GaussLine2 &gauss = cache.gauss;
     MemLine2 &mem = getMem();
 
diff --git a/tbox/src/wLine2.h b/tbox/src/wLine2.h
index 52aff488b89d9d65c2d52812255db8bdd3a5e878..400d87dac3f26eecdccf2fd9ff100c36a0e1d0cb 100644
--- a/tbox/src/wLine2.h
+++ b/tbox/src/wLine2.h
@@ -35,7 +35,7 @@ class TBOX_API Line2 : public Element
     friend class MemLine2;
     MemLine2 *pmem; ///< private memory of precalculated values
     inline MemLine2 &getMem() const;
-    inline static CacheLine2 &getCache();
+    inline static CacheLine2 &getCache(int n);
     virtual double buildJ(size_t k) const override;
     virtual void buildGradientJ(size_t k, std::vector<double> &dDetJ) const override;
     virtual std::vector<double> buildGradientV(int dim) const override;
diff --git a/tbox/src/wLine2.inl b/tbox/src/wLine2.inl
index e9e2004bee2709470661bea68ed322c5886bf4fd..704f23146100a2e1ac7115e688f82e3ee85f6205 100644
--- a/tbox/src/wLine2.inl
+++ b/tbox/src/wLine2.inl
@@ -25,8 +25,8 @@ inline MemLine2 &Line2::getMem() const
 /**
  * @brief Return cache for Line2 private methods
  */
-inline CacheLine2 &Line2::getCache()
+inline CacheLine2 &Line2::getCache(int n)
 {
-    static Lazy<CacheLine2> cache(2);
+    static Lazy<CacheLine2> cache(n);
     return cache.get();
 }
diff --git a/tbox/src/wPoint1.cpp b/tbox/src/wPoint1.cpp
index 609cc0b74b076675aa9239f9db815b001c64c009..e4cb27afa6e57dd3542f9d0d92892bd67ebd7aff 100644
--- a/tbox/src/wPoint1.cpp
+++ b/tbox/src/wPoint1.cpp
@@ -25,15 +25,13 @@ Point1::Point1(int n, Tag *_ptag, Tag *_etag,
 {
 }
 
-CachePoint1 &
-Point1::getCache()
+CachePoint1 &Point1::getCache(int n)
 {
-    static Lazy<CachePoint1> cache(2);
+    static Lazy<CachePoint1> cache(n);
     return cache.get();
 }
 
-Cache &
-Point1::getVCache() const
+Cache &Point1::getVCache() const
 {
-    return getCache();
+    return getCache(order);
 }
diff --git a/tbox/src/wPoint1.h b/tbox/src/wPoint1.h
index 2cd9517c67fe47e6fb44486d4f9940aa294671f2..fc78eb37711c4af8828be1100b83465ded5e8121 100644
--- a/tbox/src/wPoint1.h
+++ b/tbox/src/wPoint1.h
@@ -34,7 +34,7 @@ public:
     virtual ELTYPE type() const override { return ELTYPE::POINT1; }
 
 private:
-    static CachePoint1 &getCache();
+    static CachePoint1 &getCache(int n);
     virtual Cache &getVCache() const override;
 };
 
diff --git a/tbox/src/wQuad4.cpp b/tbox/src/wQuad4.cpp
index 7710081d5a143e60db459a5c206c9f3f2d9fc5fb..c44ab416db48fcde6b60c8115a704aa32c25b04e 100644
--- a/tbox/src/wQuad4.cpp
+++ b/tbox/src/wQuad4.cpp
@@ -27,7 +27,7 @@ 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().gauss.getN()));
+    pmem = new MemQuad4(*this, static_cast<int>(getCache(order).gauss.getN()));
 }
 Quad4::~Quad4()
 {
@@ -48,7 +48,7 @@ Mem &Quad4::getVMem() const
  */
 Cache &Quad4::getVCache() const
 {
-    return getCache();
+    return getCache(order);
 }
 
 /**
@@ -68,7 +68,7 @@ std::vector<Eigen::Vector3d> Quad4::computeTangents() const
  */
 std::vector<Eigen::Vector3d> Quad4::buildTangents(size_t k) const
 {
-    Eigen::MatrixXd const &dff = getCache().getDsf(k);
+    Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
 
     std::vector<Eigen::Vector3d> t(2, Eigen::Vector3d::Zero());
     size_t i = 0;
@@ -106,7 +106,7 @@ double Quad4::buildJ(size_t k) const
  */
 double Quad4::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
 {
-    Eigen::MatrixXd const &dff = getCache().getDsf(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;
@@ -124,7 +124,7 @@ double Quad4::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
  */
 Eigen::MatrixXd Quad4::buildS()
 {
-    CacheQuad4 &cache = getCache();
+    CacheQuad4 &cache = getCache(order);
     GaussQuad4 &gauss = cache.gauss;
     MemQuad4 &mem = getMem();
 
@@ -165,7 +165,7 @@ Eigen::MatrixXd Quad4::build(MATTYPE type)
  */
 double Quad4::computeV(std::vector<double> const &u, Fct0 const &fct)
 {
-    CacheQuad4 &cache = getCache();
+    CacheQuad4 &cache = getCache(order);
     GaussQuad4 &gauss = cache.gauss;
     MemQuad4 &mem = getMem();
 
diff --git a/tbox/src/wQuad4.h b/tbox/src/wQuad4.h
index 6f2e81e0e8f5b13f2039976bfcae2e4de00e4d17..fa8d1060fc08d652b1df520a6d46ce4e9c37f2f0 100644
--- a/tbox/src/wQuad4.h
+++ b/tbox/src/wQuad4.h
@@ -34,7 +34,7 @@ class TBOX_API Quad4 : public Element
 #ifndef SWIG
     friend class MemQuad4;
     MemQuad4 *pmem; ///< private memory of precalculated values
-    inline static CacheQuad4 &getCache();
+    inline static CacheQuad4 &getCache(int n);
     inline MemQuad4 &getMem() const;
     virtual double buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const override;
     virtual double buildJ(size_t k) const override;
diff --git a/tbox/src/wQuad4.inl b/tbox/src/wQuad4.inl
index 688da69d154769f7de75463118630f7a96e804af..6e23c1b10e1aadfaa77f2b3663c933037199b74f 100644
--- a/tbox/src/wQuad4.inl
+++ b/tbox/src/wQuad4.inl
@@ -25,8 +25,8 @@ inline MemQuad4 &Quad4::getMem() const
 /**
  * @brief Return cache for Quad4 private methods
  */
-inline CacheQuad4 &Quad4::getCache()
+inline CacheQuad4 &Quad4::getCache(int n)
 {
-    static Lazy<CacheQuad4> cache(2);
+    static Lazy<CacheQuad4> cache(n);
     return cache.get();
 }
diff --git a/tbox/src/wTetra4.cpp b/tbox/src/wTetra4.cpp
index 815d1a0714f07bf027d191bc8e0604618f99b1f4..fe2e1fca0fc17a338aedcf13089dd2ace39c5c0c 100644
--- a/tbox/src/wTetra4.cpp
+++ b/tbox/src/wTetra4.cpp
@@ -27,7 +27,7 @@ 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().gauss.getN()));
+    pmem = new MemTetra4(*this, static_cast<int>(getCache(order).gauss.getN()));
 }
 Tetra4::~Tetra4()
 {
@@ -48,7 +48,7 @@ Mem &Tetra4::getVMem() const
  */
 Cache &Tetra4::getVCache() const
 {
-    return getCache();
+    return getCache(order);
 }
 
 /**
@@ -58,7 +58,7 @@ Cache &Tetra4::getVCache() const
  */
 double Tetra4::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
 {
-    Eigen::MatrixXd const &dff = getCache().getDsf(k);
+    Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
 
     Eigen::Matrix3d JJ = Eigen::Matrix3d::Zero(); // temporary fixed-size matrix to efficiently compute the inverse
     size_t i = 0;
@@ -77,7 +77,7 @@ double Tetra4::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
 void Tetra4::buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eigen::MatrixXd> &dJ) const
 {
     MemTetra4 &mem = getMem();
-    Eigen::MatrixXd const &dsf = getCache().getDsf(k);
+    Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
 
     // Jacobian
     dJ.resize(4 * 3); // 4 nodes and 3 dimensions (x, y, z)
@@ -103,7 +103,7 @@ void Tetra4::buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Ei
  */
 std::vector<double> Tetra4::buildGradientV(int dim) const
 {
-    GaussTetra4 &gauss = getCache().gauss;
+    GaussTetra4 &gauss = getCache(order).gauss;
     MemTetra4 &mem = getMem();
 
     // Gradient of volume
@@ -125,7 +125,7 @@ std::vector<double> Tetra4::buildGradientV(int dim) const
  */
 Eigen::MatrixXd Tetra4::buildM()
 {
-    CacheTetra4 &cache = getCache();
+    CacheTetra4 &cache = getCache(order);
     GaussTetra4 &gauss = cache.gauss;
     MemTetra4 &mem = getMem();
 
@@ -147,7 +147,7 @@ Eigen::MatrixXd Tetra4::buildM()
  */
 Eigen::MatrixXd Tetra4::buildK(std::vector<double> const &u, Fct0 const &fct)
 {
-    CacheTetra4 &cache = getCache();
+    CacheTetra4 &cache = getCache(order);
     GaussTetra4 &gauss = cache.gauss;
     MemTetra4 &mem = getMem();
 
@@ -172,7 +172,7 @@ Eigen::MatrixXd Tetra4::buildK(std::vector<double> const &u, Fct0 const &fct)
  */
 Eigen::MatrixXd Tetra4::buildDK(std::vector<double> const &u, Fct0 const &fct)
 {
-    CacheTetra4 &cache = getCache();
+    CacheTetra4 &cache = getCache(order);
     GaussTetra4 &gauss = cache.gauss;
     MemTetra4 &mem = getMem();
 
@@ -198,7 +198,7 @@ Eigen::MatrixXd Tetra4::buildDK(std::vector<double> const &u, Fct0 const &fct)
  */
 Eigen::VectorXd Tetra4::buildDs(std::vector<double> const &u, Fct0 const &f)
 {
-    CacheTetra4 &cache = getCache();
+    CacheTetra4 &cache = getCache(order);
     GaussTetra4 &gauss = cache.gauss;
     MemTetra4 &mem = getMem();
 
@@ -216,7 +216,7 @@ Eigen::VectorXd Tetra4::buildDs(std::vector<double> const &u, Fct0 const &f)
  */
 Eigen::MatrixXd Tetra4::buildK_mechanical(Eigen::MatrixXd const &H)
 {
-    CacheTetra4 &cache = getCache();
+    CacheTetra4 &cache = getCache(order);
     GaussTetra4 &gauss = cache.gauss;
     MemTetra4 &mem = getMem();
 
@@ -254,7 +254,7 @@ Eigen::MatrixXd Tetra4::buildK_mechanical(Eigen::MatrixXd const &H)
  */
 Eigen::MatrixXd Tetra4::buildS_thermomechanical(Eigen::MatrixXd const &D)
 {
-    CacheTetra4 &cache = getCache();
+    CacheTetra4 &cache = getCache(order);
     GaussTetra4 &gauss = cache.gauss;
     MemTetra4 &mem = getMem();
 
@@ -279,7 +279,7 @@ Eigen::MatrixXd Tetra4::buildS_thermomechanical(Eigen::MatrixXd const &D)
  */
 Eigen::MatrixXd Tetra4::buildL_thermal(Eigen::MatrixXd const &C)
 {
-    CacheTetra4 &cache = getCache();
+    CacheTetra4 &cache = getCache(order);
     GaussTetra4 &gauss = cache.gauss;
     MemTetra4 &mem = getMem();
 
@@ -304,7 +304,7 @@ Eigen::MatrixXd Tetra4::buildL_thermal(Eigen::MatrixXd const &C)
  */
 void Tetra4::strain_stress(Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, Eigen::MatrixXd const &H, std::vector<double> const &u)
 {
-    CacheTetra4 &cache = getCache();
+    CacheTetra4 &cache = getCache(order);
     GaussTetra4 &gauss = cache.gauss;
     MemTetra4 &mem = getMem();
 
@@ -344,7 +344,7 @@ Eigen::VectorXd Tetra4::computeGrad(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().getDsf(k) * ue;
+    return getMem().getJinv(k) * getCache(order).getDsf(k) * ue;
 }
 
 /**
@@ -362,7 +362,7 @@ double Tetra4::computeGrad2(std::vector<double> const &u, size_t k)
  */
 double Tetra4::computeV(std::vector<double> const &u, Fct0 const &fct)
 {
-    CacheTetra4 &cache = getCache();
+    CacheTetra4 &cache = getCache(order);
     GaussTetra4 &gauss = cache.gauss;
     MemTetra4 &mem = getMem();
 
diff --git a/tbox/src/wTetra4.h b/tbox/src/wTetra4.h
index 59338935498681eb42cd58f84b0787b46903f72e..76108f233aac893885e83cccb554629d0fccd189 100644
--- a/tbox/src/wTetra4.h
+++ b/tbox/src/wTetra4.h
@@ -34,7 +34,7 @@ class TBOX_API Tetra4 : public Element
 #ifndef SWIG
     friend class MemTetra4;
     MemTetra4 *pmem; ///< private memory of precalculated values
-    inline static CacheTetra4 &getCache();
+    inline static CacheTetra4 &getCache(int n);
     inline MemTetra4 &getMem() const;
     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;
diff --git a/tbox/src/wTetra4.inl b/tbox/src/wTetra4.inl
index 4589120ea95307ee5798677e99423bd248c37b93..506df073687d813473395d1e53cb172b695fb7b6 100644
--- a/tbox/src/wTetra4.inl
+++ b/tbox/src/wTetra4.inl
@@ -25,8 +25,8 @@ inline MemTetra4 &Tetra4::getMem() const
 /**
  * @brief Return cache for Tetra4 private methods
  */
-inline CacheTetra4 &Tetra4::getCache()
+inline CacheTetra4 &Tetra4::getCache(int n)
 {
-    static Lazy<CacheTetra4> cache(2);
+    static Lazy<CacheTetra4> cache(n);
     return cache.get();
 }
diff --git a/tbox/src/wTri3.cpp b/tbox/src/wTri3.cpp
index b6493264b386dd11027807ae40ff8349194a99ed..b5a107db15fac483cf2a0e8697d036345533fa8f 100644
--- a/tbox/src/wTri3.cpp
+++ b/tbox/src/wTri3.cpp
@@ -29,7 +29,7 @@ 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().gauss.getN()));
+    pmem = new MemTri3(*this, static_cast<int>(getCache(order).gauss.getN()));
 }
 Tri3::~Tri3()
 {
@@ -50,7 +50,7 @@ Mem &Tri3::getVMem() const
  */
 Cache &Tri3::getVCache() const
 {
-    return getCache();
+    return getCache(order);
 }
 
 /**
@@ -70,7 +70,7 @@ std::vector<Eigen::Vector3d> Tri3::computeTangents() const
  */
 std::vector<Eigen::Vector3d> Tri3::buildTangents(size_t k) const
 {
-    Eigen::MatrixXd const &dff = getCache().getDsf(k);
+    Eigen::MatrixXd const &dff = getCache(order).getDsf(k);
 
     std::vector<Eigen::Vector3d> t(2, Eigen::Vector3d::Zero());
     size_t i = 0;
@@ -133,7 +133,7 @@ double Tri3::buildJ(size_t k) const
  */
 double Tri3::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
 {
-    Eigen::MatrixXd const &dff = getCache().getDsf(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;
@@ -151,7 +151,7 @@ double Tri3::buildJ(size_t k, Eigen::MatrixXd &J, Eigen::MatrixXd &iJ) const
  */
 void Tri3::buildGradientJ(size_t k, std::vector<double> &dDetJ) const
 {
-    Eigen::MatrixXd const &dsf = getCache().getDsf(k);
+    Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
 
     // Tangent vectors and determinant
     std::vector<Eigen::Vector3d> ts = this->buildTangents(k);
@@ -178,7 +178,7 @@ void Tri3::buildGradientJ(size_t k, std::vector<double> &dDetJ) const
 void Tri3::buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eigen::MatrixXd> &dJ) const
 {
     MemTri3 &mem = getMem();
-    Eigen::MatrixXd const &dsf = getCache().getDsf(k);
+    Eigen::MatrixXd const &dsf = getCache(order).getDsf(k);
 
     // Jacobian
     dJ.resize(3 * 2); // 3 nodes and 2 dimensions (x, y)
@@ -204,7 +204,7 @@ void Tri3::buildGradientJ(size_t k, std::vector<double> &dDetJ, std::vector<Eige
  */
 std::vector<double> Tri3::buildGradientV(int dim) const
 {
-    GaussTri3 &gauss = getCache().gauss;
+    GaussTri3 &gauss = getCache(order).gauss;
     MemTri3 &mem = getMem();
 
     // Gradient of volume
@@ -226,7 +226,7 @@ std::vector<double> Tri3::buildGradientV(int dim) const
  */
 Eigen::MatrixXd Tri3::buildK(std::vector<double> const &u, Fct0 const &fct)
 {
-    CacheTri3 &cache = getCache();
+    CacheTri3 &cache = getCache(order);
     GaussTri3 &gauss = cache.gauss;
     MemTri3 &mem = getMem();
 
@@ -251,7 +251,7 @@ Eigen::MatrixXd Tri3::buildK(std::vector<double> const &u, Fct0 const &fct)
  */
 Eigen::MatrixXd Tri3::buildK(std::vector<double> const &u, Fct2 const &fct, bool fake)
 {
-    CacheTri3 &cache = getCache();
+    CacheTri3 &cache = getCache(order);
     GaussTri3 &gauss = cache.gauss;
 
     Eigen::MatrixXd K = Eigen::Matrix3d::Zero();
@@ -289,7 +289,7 @@ Eigen::MatrixXd Tri3::buildK(std::vector<double> const &u, Fct2 const &fct, bool
  */
 Eigen::MatrixXd Tri3::buildK_mechanical(Eigen::MatrixXd const &H)
 {
-    CacheTri3 &cache = getCache();
+    CacheTri3 &cache = getCache(order);
     GaussTri3 &gauss = cache.gauss;
     MemTri3 &mem = getMem();
 
@@ -325,7 +325,7 @@ Eigen::MatrixXd Tri3::buildK_mechanical(Eigen::MatrixXd const &H)
  */
 Eigen::MatrixXd Tri3::buildDK(std::vector<double> const &u, Fct0 const &fct)
 {
-    CacheTri3 &cache = getCache();
+    CacheTri3 &cache = getCache(order);
     GaussTri3 &gauss = cache.gauss;
     MemTri3 &mem = getMem();
 
@@ -351,7 +351,7 @@ Eigen::MatrixXd Tri3::buildDK(std::vector<double> const &u, Fct0 const &fct)
  */
 Eigen::VectorXd Tri3::builds(std::vector<double> const &u, Fct0 const &f)
 {
-    CacheTri3 &cache = getCache();
+    CacheTri3 &cache = getCache(order);
     GaussTri3 &gauss = cache.gauss;
     MemTri3 &mem = getMem();
 
@@ -369,7 +369,7 @@ Eigen::VectorXd Tri3::builds(std::vector<double> const &u, Fct0 const &f)
  */
 Eigen::VectorXd Tri3::builds(std::vector<double> const &u, Fct1 const &f)
 {
-    CacheTri3 &cache = getCache();
+    CacheTri3 &cache = getCache(order);
     GaussTri3 &gauss = cache.gauss;
     MemTri3 &mem = getMem();
 
@@ -387,7 +387,7 @@ Eigen::VectorXd Tri3::builds(std::vector<double> const &u, Fct1 const &f)
  */
 Eigen::VectorXd Tri3::buildDs(std::vector<double> const &u, Fct0 const &f)
 {
-    CacheTri3 &cache = getCache();
+    CacheTri3 &cache = getCache(order);
     GaussTri3 &gauss = cache.gauss;
     MemTri3 &mem = getMem();
 
@@ -405,7 +405,7 @@ Eigen::VectorXd Tri3::buildDs(std::vector<double> const &u, Fct0 const &f)
  */
 Eigen::MatrixXd Tri3::buildS()
 {
-    CacheTri3 &cache = getCache();
+    CacheTri3 &cache = getCache(order);
     GaussTri3 &gauss = cache.gauss;
     MemTri3 &mem = getMem();
 
@@ -432,7 +432,7 @@ Eigen::VectorXd Tri3::computeGrad(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().getDsf(k) * ue;
+    return getMem().getJinv(k) * getCache(order).getDsf(k) * ue;
 }
 
 /**
@@ -450,7 +450,7 @@ double Tri3::computeGrad2(std::vector<double> const &u, size_t k)
  */
 Eigen::VectorXd Tri3::computeqint(std::vector<double> const &u, Fct2 const &fct)
 {
-    CacheTri3 &cache = getCache();
+    CacheTri3 &cache = getCache(order);
     GaussTri3 &gauss = cache.gauss;
     MemTri3 &mem = getMem();
 
@@ -476,7 +476,7 @@ Eigen::VectorXd Tri3::computeqint(std::vector<double> const &u, Fct2 const &fct)
  */
 Eigen::VectorXd Tri3::computeqV(std::vector<double> const &u, Fct2 const &fct)
 {
-    CacheTri3 &cache = getCache();
+    CacheTri3 &cache = getCache(order);
     GaussTri3 &gauss = cache.gauss;
     MemTri3 &mem = getMem();
 
@@ -501,7 +501,7 @@ Eigen::VectorXd Tri3::computeqV(std::vector<double> const &u, Fct2 const &fct)
  */
 double Tri3::computeV(std::vector<double> const &u, Fct0 const &fct)
 {
-    CacheTri3 &cache = getCache();
+    CacheTri3 &cache = getCache(order);
     GaussTri3 &gauss = cache.gauss;
     MemTri3 &mem = getMem();
 
@@ -519,12 +519,12 @@ double Tri3::computeV(std::vector<double> const &u, Fct0 const &fct)
  */
 Eigen::MatrixXd Tri3::computeV(std::vector<double> const &u, Fct2 const &fct)
 {
-    CacheTri3 &cache = getCache();
+    CacheTri3 &cache = getCache(order);
     GaussTri3 &gauss = cache.gauss;
     MemTri3 &mem = getMem();
 
     // Gauss integration
-    Eigen::MatrixXd out = Eigen::MatrixXd::Zero(2, 2); /* @todo out should be resized to match fk size */
+    Eigen::MatrixXd out = Eigen::Matrix2d::Zero(); /* @todo out should be resized to match fk size */
     for (size_t k = 0; k < gauss.getN(); ++k)
     {
         // Function evaluation
diff --git a/tbox/src/wTri3.h b/tbox/src/wTri3.h
index 3fe60c3c690f60dc393be00557258218e0f77177..7255a33cd1bc7cf21fca3329cf5e71d359ac4de9 100644
--- a/tbox/src/wTri3.h
+++ b/tbox/src/wTri3.h
@@ -35,7 +35,7 @@ class TBOX_API Tri3 : public Element
     friend class MemTri3;
     MemTri3 *pmem; ///< private memory of precalculated values
     inline MemTri3 &getMem() const;
-    inline static CacheTri3 &getCache();
+    inline static CacheTri3 &getCache(int order);
     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;
diff --git a/tbox/src/wTri3.inl b/tbox/src/wTri3.inl
index 57f2a705cd8d737462d9bdc5a2209fe1f8a97b00..23df4f1e2030610350cc0d48008591d08540dd2f 100644
--- a/tbox/src/wTri3.inl
+++ b/tbox/src/wTri3.inl
@@ -25,8 +25,8 @@ inline MemTri3 &Tri3::getMem() const
 /**
  * @brief Return cache for Tri3 private methods
  */
-inline CacheTri3 &Tri3::getCache()
+inline CacheTri3 &Tri3::getCache(int n)
 {
-    static Lazy<CacheTri3> cache(2);
+    static Lazy<CacheTri3> cache(n);
     return cache.get();
 }