Update waves

Created by: acrovato

Context

This pull request brings the last developments from the flow module as well as some upgrade to waves.

Major changes in flow

• flow is now able to compute 3D flows
• The Jacobian (tangent) matrix for the Netwton method is now fully analytical
• The adjacent elements identification is now based on ```std::unordered_set````

Major changes in waves

• tbox and fwk are now included in a namespace
• waves is now fully included in a namespace
• All elements now have a private memory (Mem) which can be publicly accessed
• All elements now have a cache including Gauss points (Gauss) which can be publicly accessed
• MshDeformation now uses a linear elasticity solver
• MshCrack class implemented (creates duplicate internal boundaries in meshes)

Test

The ctest battery has been run on gaston, 100% tests passed under 110s.

Future work

• Use automated indentation rules for the whole code
• Move methods using VTK outside of MshData
• Clean tbox (mainly MshData, Tri3, Tera4 and Hex8)

flow/src/wKuttaElement.cpp

@@ -3,6 +3,7 @@
 #include "wGaussLine2.h"
 #include "wSfLine2.h"
 #include "wSfTri3.h"
+#include "wMem.h"
 using namespace flow;
 
 KuttaElement::KuttaElement(int _no, Element * &_surUpE, Element * &_surLwE,
@@ -25,25 +26,6 @@ KuttaElement::KuttaElement(int _no, Element * &_surUpE, Element * &_surLwE,
     nColLw = dffVolLw[0].size();
 }
 
-/**
- * @brief Update the Jacobian
- * @authors Adrien Crovato
- */
-void KuttaElement::update()
-{
-    // Build determinant of Jacobian (surface)
-    dtmJSur.resize(nGP);
-    for (int k = 0; k < nGP; ++k)
-        dtmJSur[k] = buildSurJ(k, surUpE);
-    // Build inverse of Jacobian (volume)
-    invJVolUp.resize(nGP);
-    invJVolLw.resize(nGP);
-    for (int k = 0; k < nGP; ++k) {
-        buildVolJ(k, invJVolUp[k], volUpE, dffVolUp);
-        buildVolJ(k, invJVolLw[k], volLwE, dffVolLw);
-    }
-}
-
 /**
  * @brief Compute stiffness matrices associated to Kutta
  * @authors Adrien Crovato
@@ -77,16 +59,18 @@ void KuttaElement::buildK(gmm::dense_matrix<double> &Kup, gmm::dense_matrix<doub
 
     // TODO: To be displaced outside
     int NDIM = 2;  // dimension of volume element
+    Mem &surMem = surUpE->getVMem();
+    Mem &volUpMem = volUpE->getVMem();
+    Mem &volLwMem = volLwE->getVMem();
     gmm::clear(Kup);
     gmm::clear(Klw);
-    update();
 
     // Build
     for (int k = 0; k < nGP; ++k) {
 
         // fill N'
         double w = wSur[k];
-        double dtm = dtmJSur[k];
+        double dtm = surMem.getDetJ(k);
         gmm::dense_matrix<double> N(nRow,1);
         for (int i = 0; i < nRow; ++i) {
             N(i,0) = ffSur[k][i];
@@ -121,8 +105,8 @@ void KuttaElement::buildK(gmm::dense_matrix<double> &Kup, gmm::dense_matrix<doub
 
         // A = inv(J)*dN
         gmm::dense_matrix<double> Aup(NDIM,nColUp), Alw(NDIM,nColLw);
-        gmm::mult(invJVolUp[k], dNup, Aup);
-        gmm::mult(invJVolLw[k], dNlw, Alw);
+        gmm::mult(volUpMem.getJinv(k), dNup, Aup);
+        gmm::mult(volLwMem.getJinv(k), dNlw, Alw);
 
         // B = V*[inv(J)*dN]
         gmm::dense_matrix<double> Bup(1,nColUp), Blw(1,nColLw);
@@ -165,15 +149,16 @@ void KuttaElement::buildNK(gmm::dense_matrix<double> &Kup, gmm::dense_matrix<dou
  */
 void KuttaElement::buildNs(std::vector<double> &bUp, std::vector<double> &bLw, std::vector<double> const &u)
 {
+    Mem &surMem = surUpE->getVMem();
     gmm::clear(bUp);
     gmm::clear(bLw);
-    update();
+
     // Build
     for (int k = 0; k < nGP; ++k) {
 
         // fill N'
         double w = wSur[k];
-        double dtm = dtmJSur[k];
+        double dtm = surMem.getDetJ(k);
         gmm::dense_matrix<double> N(nRow,1);
         for (int i = 0; i < nRow; ++i) {
             N(i,0) = ffSur[k][i];
@@ -227,32 +212,21 @@ void KuttaElement::buildSurF(Element * const &elem)
     }
     // dff = [vA vB]
     // v = [dxiN'1 dxiN'2]
-    dffSur.resize(gauss.p.size());
-    for (size_t k = 0; k < gauss.p.size(); ++k)
-    {
-        Pt ksi(gauss.p[k]);
-        dffSur[k].resize(2);
-        for (int i = 0; i < 2; ++i)
-        {
-            dffSur[k][i].resize(1);
-            sf.evalD(i, ksi, dffSur[k][i]);
-        }
-    }
+    //dffSur.resize(gauss.p.size());
+    //for (size_t k = 0; k < gauss.p.size(); ++k)
+    //{
+    //    Pt ksi(gauss.p[k]);
+    //    dffSur[k].resize(2);
+    //    for (int i = 0; i < 2; ++i)
+    //    {
+    //        dffSur[k][i].resize(1);
+    //        sf.evalD(i, ksi, dffSur[k][i]);
+    //    }
+    //}
     wSur.resize(nGP);
     wSur = gauss.w;
 }
 
-/**
- * @brief Compute surface Jacobian
- * @authors Adrien Crovato
- */
-double KuttaElement::buildSurJ(int k, Element * const &elem)
-{
-    if (elem->type() != ELTYPE::LINE2) throw std::runtime_error("KuttaElement not implemented\n");
-
-    return elem->buildJ(k);
-}
-
 /**
  * @brief Compute volume shape functions and derivatives
  * @authors Adrien Crovato
@@ -282,20 +256,6 @@ void KuttaElement::buildVolF(Element * const &elem, std::vector<std::vector<std:
     }
 }
 
-/**
- * @brief Compute inverse of volume Jacobian
- * @authors Adrien Crovato
- */
-void KuttaElement::buildVolJ(int k, gmm::dense_matrix<double> &J, Element * const &elem, std::vector<std::vector<std::vector<double>>> const &dffVol)
-{
-    if (elem->type() != ELTYPE::TRI3) throw std::runtime_error("KuttaElement not implemented\n");
-
-    J.resize(2, 2);
-    elem->buildJ(k, J);
-    // J = inv(J)
-    gmm::lu_inverse(J);
-}
-
 void KuttaElement::write(std::ostream &out) const
 {
     out << "KuttaElement #" << no