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/wSolver.cpp

@@ -94,6 +94,8 @@ Solver::Solver(std::shared_ptr<Problem> _pbl) : pbl(_pbl)
         iIC->apply(phi);
     for(auto dBC : pbl->dBCs)
         dBC->apply(phi);
+    // Find initial best upwind element
+    findUpwd();
 }
 Solver::~Solver()
 {
@@ -121,6 +123,7 @@ bool Solver::run()
     // Initialize solver loop
     int itC = 0; // iteration counter
     int avC = 0; // adaptive viscosity counter
+    bool rUpwd = true; // re-find best upwind element
 
     // Compute initial residual
     buildRes(rPhi);
@@ -151,8 +154,10 @@ bool Solver::run()
     do
     {
         // Viscosity damping
-        if (relRes < avThrsh)
+        if (relRes < avThrsh && avC < 2) {
             avC++;
+            rUpwd = true;
+        }
         switch(avC) {
             case 0 : {  mCOv = 0.92;
                         muCv = 2.;
@@ -196,6 +201,11 @@ bool Solver::run()
         {}
         absRes = gmm::vect_norm2(rPhi);
         relRes = absRes / resInit;
+        // Re-find best upwind element
+        if (rUpwd) {
+            findUpwd();
+            rUpwd = false;
+        }
 
         // Compute aerodynamic coefficient on boundaries
         Cl = 0.; Cd = 0.;
@@ -266,8 +276,7 @@ void Solver::buildJac(gmm::csr_matrix<double> &J)
         // Current element
         Element * e = p.first;
         //Upwind element
-        size_t idx = findBestUpwd(p);
-        Element * eU = p.second[idx];
+        Element * eU = p.second[0];
 
         // Merge nodes of current and upwind elements (and remove duplicate)
         std::vector<Node *> nodeU;
@@ -405,8 +414,7 @@ void Solver::buildRes(std::vector<double> &R)
         // Current element
         Element * e = p.first;
         // Upwind element
-        size_t idx = findBestUpwd(p);
-        Element * eU = p.second[idx];
+        Element * eU = p.second[0];
 
         // Switching function and bias
         double mach = fluid->mach.eval(*e, phi, 0);
@@ -549,31 +557,37 @@ void Solver::compute()
  * @brief Find upwind element which is best aligned with current velocity vector
  * @authors Adrien Crovato
  */
-size_t Solver::findBestUpwd(std::pair<Element *, std::vector<Element *>> &p)
+void Solver::findUpwd()
 {
-    // Velocity
-    std::vector<double> gradPhi(pbl->nDim);
-    p.first->computeGrad(phi, gradPhi, 0);
-    Pt vel_;
-    for (size_t i = 0; i < gradPhi.size(); ++i)
-        vel_.x[i] = -gradPhi[i];
-    vel_.normalize();
-    // Best alignement strategy
-    // -> find (cgU-cgE) \dot vel_ closest to one
-    double SPbest = -1;
-    size_t idx = 0;
-    for (size_t i = 0; i < p.second.size(); ++i) {
-        Pt cgU = (p.second[i])->getVMem().getCg();
-        Pt cgE = (p.first)->getVMem().getCg();
-        double SP = (cgU - cgE).normalized() * vel_;
-        if (fabs(SP-1) < fabs(SPbest-1)) {
-            idx = i;
-            SPbest = SP;
+    tbb::parallel_do(pbl->medium->adjMap.begin(), pbl->medium->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> &p)
+	{
+        // Velocity
+        std::vector<double> gradPhi(pbl->nDim);
+        p.first->computeGrad(phi, gradPhi, 0);
+        Pt vel_;
+        for (size_t i = 0; i < gradPhi.size(); ++i)
+            vel_.x[i] = -gradPhi[i];
+        vel_.normalize();
+        // Best alignement strategy
+        // -> find (cgU-cgE) \dot vel_ closest to one
+        double SPbest = -1;
+        size_t idx = 0;
+        for (size_t i = 0; i < p.second.size(); ++i) {
+            Pt cgU = (p.second[i])->getVMem().getCg();
+            Pt cgE = (p.first)->getVMem().getCg();
+            double SP = (cgU - cgE).normalized() * vel_;
+            if (fabs(SP-1) < fabs(SPbest-1)) {
+                idx = i;
+                SPbest = SP;
+            }
+            else
+                continue;
         }
-        else
-            continue;
-    }
-    return idx;
+        // Swap (best upwind element is moved at first emplacement)
+        Element * tmp = p.second[0];
+		p.second[0] = p.second[idx];
+		p.second[idx] = tmp;
+    });
 }
 
 /**