From c781978162d4d8a7967fe18764bc7c0eb502a396 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Fri, 4 Jun 2021 14:44:18 +0200
Subject: [PATCH] Add variable controling order of quadrature. Add quadrature
 rule for Tetra4 of degree 1.

---
 tbox/src/wCache.h         |  2 +-
 tbox/src/wCacheHex8.h     |  2 +-
 tbox/src/wCacheLine2.h    |  2 +-
 tbox/src/wCachePoint1.h   |  2 +-
 tbox/src/wCacheQuad4.h    |  2 +-
 tbox/src/wCacheTetra4.h   |  2 +-
 tbox/src/wCacheTri3.h     |  2 +-
 tbox/src/wElement.cpp     |  3 ++-
 tbox/src/wElement.h       |  2 ++
 tbox/src/wGaussTetra4.cpp |  9 +++++++++
 tbox/src/wHex8.cpp        | 20 +++++++++----------
 tbox/src/wHex8.h          |  2 +-
 tbox/src/wHex8.inl        |  4 ++--
 tbox/src/wLine2.cpp       | 16 +++++++--------
 tbox/src/wLine2.h         |  2 +-
 tbox/src/wLine2.inl       |  4 ++--
 tbox/src/wPoint1.cpp      | 10 ++++------
 tbox/src/wPoint1.h        |  2 +-
 tbox/src/wQuad4.cpp       | 12 +++++------
 tbox/src/wQuad4.h         |  2 +-
 tbox/src/wQuad4.inl       |  4 ++--
 tbox/src/wTetra4.cpp      | 30 ++++++++++++++--------------
 tbox/src/wTetra4.h        |  2 +-
 tbox/src/wTetra4.inl      |  4 ++--
 tbox/src/wTri3.cpp        | 42 +++++++++++++++++++--------------------
 tbox/src/wTri3.h          |  2 +-
 tbox/src/wTri3.inl        |  4 ++--
 27 files changed, 100 insertions(+), 90 deletions(-)

diff --git a/tbox/src/wCache.h b/tbox/src/wCache.h
index e5dab88f..f20921ce 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 b0f6e9a7..212dba30 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 1d0f160f..a912013a 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 d14baf72..9d5f290b 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 11cc0210..49a2fb1e 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 b56b1b1b..4d1062ff 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 da73a5d7..c7491eb7 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 c3cd1127..a1b3d561 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 250884ef..48d60c92 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 5bdff84d..0e1a4e35 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 792c5647..da1332f5 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 7850cc98..1c563df0 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 aa6e937f..6912e483 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 00dff4ad..986904ef 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 52aff488..400d87da 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 e9e2004b..704f2314 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 609cc0b7..e4cb27af 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 2cd9517c..fc78eb37 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 7710081d..c44ab416 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 6f2e81e0..fa8d1060 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 688da69d..6e23c1b1 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 815d1a07..fe2e1fca 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 59338935..76108f23 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 4589120e..506df073 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 b6493264..b5a107db 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 3fe60c3c..7255a33c 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 57f2a705..23df4f1e 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();
 }
-- 
GitLab