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

@@ -268,52 +268,73 @@ void Solver::buildJac(gmm::csr_matrix<double> &J)
         //Upwind element
         Element * eU = p.second[0];
 
-        // Merge nodes of current and upwind elements (and remove duplicate)
-        std::vector<Node *> nodeU;
-        for (auto n : e->nodes)
-            nodeU.push_back(n);
-        for (auto n : eU->nodes)
-            nodeU.push_back(n);
-        std::sort (nodeU.begin(), nodeU.end());
-        auto iterator = std::unique(nodeU.begin(), nodeU.end());
-        nodeU.resize(std::distance(nodeU.begin(), iterator));
-
-        // Perturbed solution
-        std::vector<double> phiPlu = phi;
-        std::vector<double> phiMin = phi;
-        double phiTmp;
-        for (size_t jj = 0; jj < nodeU.size(); ++jj) {
-            Node * nodj = nodeU[jj];
-            // Perturb the solution
-            phiTmp = phi[nodj->row];
-            phiPlu[nodj->row] += epsNR;
-            phiMin[nodj->row] -= epsNR;
-
-            // Elementary flux vector
-            std::vector<double> bePlu(e->nodes.size());
-            std::vector<double> beMin(e->nodes.size());
-
-            // Switching function and biais
-            double machPlu = fluid->mach.eval(*e, phiPlu, 0);
-            double machMin = fluid->mach.eval(*e, phiMin, 0);
-            double muPlu = std::max(muCv*(1-(mCOv*mCOv)/(machPlu*machPlu)) ,0.);
-            double muMin = std::max(muCv*(1-(mCOv*mCOv)/(machMin*machMin)), 0.);
-            double biasPlu = muPlu * (fluid->rho.eval(*e, phiPlu, 0) - fluid->rho.eval(*eU, phiPlu, 0));
-            double biasMin = muMin * (fluid->rho.eval(*e, phiMin, 0) - fluid->rho.eval(*eU, phiMin, 0));
-            e->buildNs(bePlu, phiPlu, fluid->rho, biasPlu);
-            e->buildNs(beMin, phiMin, fluid->rho, biasMin);
-
+        // If flow is:
+        // - subsonic, no upwinding, analytical Jacobian
+        // - supersonic, add upwinding, numerical jacobian
+        double mach = fluid->mach.eval(*e, phi, 0);
+        if (mach < mCOv) {
+            gmm::dense_matrix<double> Ae(e->nodes.size(), e->nodes.size());
+            e->buildNK(Ae, phi, fluid->rho);
             // Assembly
-            tbb::spin_mutex::scoped_lock lock(mutex);
-            for(size_t ii = 0; ii < e->nodes.size(); ++ii)
-            {
-                Node *nodi = e->nodes[ii];
-                A(nodi->row, nodj->row) += (bePlu[ii] - beMin[ii]) / (2*epsNR);
+	        tbb::spin_mutex::scoped_lock lock(mutex);
+	        for(size_t ii=0; ii<e->nodes.size(); ++ii)
+	        {
+	        	Node *nodi = e->nodes[ii];
+	        	for(size_t jj=0; jj<e->nodes.size(); ++jj)
+	        	{
+	        		Node *nodj = e->nodes[jj];
+	        		A(nodi->row,nodj->row) += Ae(ii,jj);
+	        	}
+	        }
+        }
+        else {
+            // Merge nodes of current and upwind elements (and remove duplicate)
+            std::vector<Node *> nodeU;
+            for (auto n : e->nodes)
+                nodeU.push_back(n);
+            for (auto n : eU->nodes)
+                nodeU.push_back(n);
+            std::sort (nodeU.begin(), nodeU.end());
+            auto iterator = std::unique(nodeU.begin(), nodeU.end());
+            nodeU.resize(std::distance(nodeU.begin(), iterator));
+
+            // Perturbed solution
+            std::vector<double> phiPlu = phi;
+            std::vector<double> phiMin = phi;
+            double phiTmp;
+            for (size_t jj = 0; jj < nodeU.size(); ++jj) {
+                Node * nodj = nodeU[jj];
+                // Perturb the solution
+                phiTmp = phi[nodj->row];
+                phiPlu[nodj->row] += epsNR;
+                phiMin[nodj->row] -= epsNR;
+
+                // Elementary flux vector
+                std::vector<double> bePlu(e->nodes.size());
+                std::vector<double> beMin(e->nodes.size());
+
+                // Switching function and biais
+                double machPlu = fluid->mach.eval(*e, phiPlu, 0);
+                double machMin = fluid->mach.eval(*e, phiMin, 0);
+                double muPlu = std::max(muCv*(1-(mCOv*mCOv)/(machPlu*machPlu)) ,0.);
+                double muMin = std::max(muCv*(1-(mCOv*mCOv)/(machMin*machMin)), 0.);
+                double biasPlu = muPlu * (fluid->rho.eval(*e, phiPlu, 0) - fluid->rho.eval(*eU, phiPlu, 0));
+                double biasMin = muMin * (fluid->rho.eval(*e, phiMin, 0) - fluid->rho.eval(*eU, phiMin, 0));
+                e->buildNs(bePlu, phiPlu, fluid->rho, biasPlu);
+                e->buildNs(beMin, phiMin, fluid->rho, biasMin);
+
+                // Assembly
+                tbb::spin_mutex::scoped_lock lock(mutex);
+                for(size_t ii = 0; ii < e->nodes.size(); ++ii)
+                {
+                    Node *nodi = e->nodes[ii];
+                    A(nodi->row, nodj->row) += (bePlu[ii] - beMin[ii]) / (2*epsNR);
+                }
+
+                // Reset the solution
+                phiPlu[nodj->row] = phiTmp;
+                phiMin[nodj->row] = phiTmp;
             }
-
-            // Reset the solution
-            phiPlu[nodj->row] = phiTmp;
-            phiMin[nodj->row] = phiTmp;
         }
 	});