dart/src/wAssign.h

@@ -83,4 +83,4 @@ public:
 
 } // namespace dart
 
-#endif //WASSIGN_H
+#endif // WASSIGN_H

dart/src/wBlowing.cpp

@@ -59,14 +59,6 @@ void Blowing::create()
     }
 }
 
-/**
- * @brief Set blowing velocity
- */
-void Blowing::setU(size_t i, double ue)
-{
-    uE[i].second->update(ue);
-}
-
 void Blowing::write(std::ostream &out) const
 {
     out << " Blowing boundary condition on " << *tag << std::endl;

dart/src/wBlowing.h

@@ -41,7 +41,8 @@ public:
     Blowing(std::shared_ptr<tbox::MshData> _msh, std::string const &names);
     virtual ~Blowing();
 
-    void setU(size_t i, double ue);
+    void setU(size_t i, double ue) { uE[i].second->update(ue); }
+    double getU(size_t i) const { return uE[i].second->eval(*uE[i].first, std::vector<double>(), 0); }
 
 #ifndef SWIG
     virtual void write(std::ostream &out) const override;
@@ -53,4 +54,4 @@ private:
 
 } // namespace dart
 
-#endif //WBLOWING_H
+#endif // WBLOWING_H

dart/src/wBlowingResidual.cpp

@@ -27,7 +27,7 @@ using namespace dart;
 /**
  * @brief Build the residual vector, on one boundary element
  *
- * b = \int psi * vn dV
+ * b = \int psi * vn dS
  */
 Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> const &phi, F0El const &f)
 {
@@ -44,3 +44,44 @@ Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> con
         b += cache.getSf(k) * vn * gauss.getW(k) * e.getDetJ(k);
     return b;
 }
+
+/**
+ * @brief Build the gradient of the residual vector with respect to the blowing velocity, on one blowing boundary element
+ *
+ * b = \int psi dS
+ */
+Eigen::VectorXd BlowingResidual::buildGradientBlowing(tbox::Element const &e)
+{
+    // Get pre-computed values
+    Cache &cache = e.getVCache();
+    Gauss &gauss = cache.getVGauss();
+
+    // Build velocity gradient
+    Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
+    for (size_t k = 0; k < gauss.getN(); ++k)
+        b += cache.getSf(k) * gauss.getW(k) * e.getDetJ(k);
+    return b;
+}
+
+/**
+ * @brief Build the gradient of the residual vector with respect to the mesh, on one blowing boundary element
+ *
+ * A = \int psi * vn ddS
+ */
+Eigen::MatrixXd BlowingResidual::buildGradientMesh(tbox::Element const &e, std::vector<double> const &phi, F0El const &f, int const nDim)
+{
+    // Get pre-computed values
+    Cache &cache = e.getVCache();
+    Gauss &gauss = cache.getVGauss();
+
+    // Build velocity gradient
+    double vn = f.eval(e, phi, 0);
+    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(e.nodes.size(), nDim * e.nodes.size());
+    for (size_t k = 0; k < gauss.getN(); ++k)
+    {
+        std::vector<double> dDetJ = e.getGradDetJ(k);
+        for (size_t i = 0; i < dDetJ.size(); i++)
+            A.col(i) += cache.getSf(k) * vn * gauss.getW(k) * dDetJ[i];
+    }
+    return A;
+}

dart/src/wBlowingResidual.h

@@ -33,7 +33,9 @@ class DART_API BlowingResidual
 public:
     // Newton
     static Eigen::VectorXd build(tbox::Element const &e, std::vector<double> const &phi, F0El const &f);
+    static Eigen::VectorXd buildGradientBlowing(tbox::Element const &e);
+    static Eigen::MatrixXd buildGradientMesh(tbox::Element const &e, std::vector<double> const &phi, F0El const &f, int const nDim);
 };
 
 } // namespace dart
-#endif //WBLOWINGRESIDUAL_H
+#endif // WBLOWINGRESIDUAL_H

dart/src/wBody.cpp

@@ -22,8 +22,6 @@
 #include "wElement.h"
 #include "wNode.h"
 #include <unordered_set>
-#include <iomanip>
-#include <fstream>
 using namespace tbox;
 using namespace dart;
 
@@ -111,84 +109,7 @@ void Body::create()
     nLoads.resize(nodes.size(), Eigen::Vector3d::Zero());
 }
 
-/**
- * @brief Save nodal data for further post-processing
- */
-void Body::save(std::string const &name, Results const &res)
-{
-    // Write to file
-    std::cout << "writing file: " << name + ".dat"
-              << "... " << std::flush;
-    std::ofstream outfile;
-    outfile.open(name + ".dat");
-    // Header
-    outfile << "$Body data" << std::endl;
-    // Aerodynamic coefficients
-    outfile << "$Aerodynamic coefficients" << std::endl;
-    outfile << "!" << std::fixed
-            << std::setw(14) << "Cl"
-            << std::setw(15) << "Cd"
-            << std::setw(15) << "Cs"
-            << std::setw(15) << "Cm"
-            << std::endl;
-    outfile << std::fixed
-            << std::setw(15) << Cl
-            << std::setw(15) << Cd
-            << std::setw(15) << Cs
-            << std::setw(15) << Cm
-            << std::endl;
-    // Elements (connectvity)
-    outfile << "$Elements" << std::endl;
-    outfile << groups[0]->tag->elems.size() << std::endl;
-    for (auto e : groups[0]->tag->elems)
-    {
-        outfile << std::fixed
-                << std::setw(10) << e->no;
-        for (auto n : e->nodes)
-            outfile << std::setw(10) << n->no;
-        outfile << std::endl;
-    }
-    //Nodes (data)
-    outfile << "$Nodal data" << std::endl;
-    outfile << nodes.size() << std::endl;
-    outfile << "!" << std::fixed
-            << std::setw(9) << "no"
-            << std::setw(15) << "x"
-            << std::setw(15) << "y"
-            << std::setw(15) << "z";
-    for (auto &p : res.scalars_at_nodes)
-        outfile << std::setw(15) << p.first;
-    for (auto &p : res.vectors_at_nodes)
-    {
-        outfile << std::setw(14) << p.first << "x"
-                << std::setw(14) << p.first << "y"
-                << std::setw(14) << p.first << "z";
-    }
-    outfile << std::endl;
-    for (auto n : nodes)
-    {
-        outfile << std::fixed
-                << std::setw(10) << n->no
-                << std::scientific << std::setprecision(6)
-                << std::setw(15) << n->pos(0)
-                << std::setw(15) << n->pos(1)
-                << std::setw(15) << n->pos(2);
-        for (auto &p : res.scalars_at_nodes)
-            outfile << std::setw(15) << (*(p.second))[n->row];
-        for (auto &p : res.vectors_at_nodes)
-            outfile << std::setw(15) << (*(p.second))[n->row](0)
-                    << std::setw(15) << (*(p.second))[n->row](1)
-                    << std::setw(15) << (*(p.second))[n->row](2);
-        outfile << std::endl;
-    }
-    // Footer
-    outfile << std::endl;
-    // Close file
-    outfile.close();
-    std::cout << "done" << std::endl;
-}
-
 void Body::write(std::ostream &out) const
 {
     out << *groups[0]->tag << " is a Body immersed in " << *groups[1]->tag << std::endl;
-}
\ No newline at end of file
+}

dart/src/wBody.h

@@ -49,8 +49,6 @@ public:
     Body(std::shared_ptr<tbox::MshData> _msh, std::vector<std::string> const &names);
     virtual ~Body() { std::cout << "~Body()\n"; }
 
-    void save(std::string const &name, tbox::Results const &res);
-
 #ifndef SWIG
     virtual void write(std::ostream &out) const override;
 #endif
@@ -61,4 +59,4 @@ private:
 
 } // namespace dart
 
-#endif //WBODY_H
+#endif // WBODY_H

dart/src/wF0El.cpp

@@ -28,153 +28,249 @@ void F0ElC::update(double _v)
 /**
  * @brief Evaluate the constant
  */
-double F0ElC::eval(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElC::eval(Element const &e, std::vector<double> const &phi, size_t k) const
 {
     return v;
 }
 /**
  * @brief Evaluate the gradient of the constant
  */
-double F0ElC::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElC::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
 {
     return 0;
 }
 
 /**
- * @brief Evaluate the nonlinear density (constant over element)
+ * @brief Precompute values for Padé approximation
+ *
+ * The input _mC2 is the square of the critical Mach number (used to compute
+ * the critical velocity).
  */
-double F0ElRho::eval(Element const &e, std::vector<double> const &u, size_t k) const
+F0ElSpeedSound::F0ElSpeedSound(double _mInf, double _mC2) : F0El(), gamma(1.4), mInf(_mInf)
 {
-    double gradU = e.computeGradient(u, k).norm(); // norm of velocity
-
-    if (gradU < gradUC) // true density
-    {
-        // particularized: pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), 1 / (gamma - 1));
-        double a = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU * gradU);
-        return sqrt(a * a * a * a * a);
-    }
+    uC = sqrt(_mC2 * (1 / (mInf * mInf) + (gamma - 1) / 2) / (1 + (gamma - 1) / 2 * _mC2)); // critical velocity
+    aC = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - uC * uC);                               // a^2(uC)
+    beta = 0 * ((gamma - 1) * mInf * mInf * uC) / aC;                                       // da^2/du(uC) / -a^2(uC)
+}
+/**
+ * @brief Evaluate the base of the isentropic density
+ *
+ * a^2 = [1 + (gamma-1)/2*Mi^2*(1-u^2)], u <= uC
+ * a^2 = a^2(uC) / (1-(da^2/du(uC)/a^2(uC))*(u/uC-1)), u > uC
+ */
+double F0ElSpeedSound::eval(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+    double u = e.computeGradient(phi, k).norm(); // norm of velocity
+    if (u <= uC)
+        return 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u * u);
     else
-    {
-        double a = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradUC * gradUC);
-        double beta = mInf * mInf * gradUC / a;
-        return sqrt(a * a * a * a * a) / (1 + beta * (gradU / gradUC - 1)); // limited density
-    }
+        return aC / (1 + beta * (u / uC - 1));
 }
 /**
- * @brief Evaluate the nonlinear density derivative (constant over element)
+ * @brief Evaluate the gradient of the base of the isentropic density
+ *
+ * da^2 = -(gamma-1)*Mi^2 * u, u <= uC
+ * da^2 = -a^2(uC) / (1-(da^2/du(uC)/a^2(uC))*(u/uC-1))^2 * (da^2/du(uC)/a^2(uC)) / uC, u > uC
+ * Computed value has to be multiplied by u to recover actual gradient
  */
-double F0ElRho::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElSpeedSound::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
 {
-    double gradU = e.computeGradient(u, k).norm(); // norm of velocity
+    double u = e.computeGradient(phi, k).norm(); // norm of velocity
+    if (u <= uC)
+        return -(gamma - 1) * mInf * mInf;
+    else
+        return -aC / ((1 + beta * (u / uC - 1)) * (1 + beta * (u / uC - 1))) * beta / uC / u;
+}
 
-    if (gradU < gradUC) // true density
-    {
-        //particularized: -mInf * mInf * pow(rho, 2 - gamma);
-        double a = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU * gradU);
-        return -mInf * mInf * sqrt(a * a * a); // has to be multiplied by grad u to recover true derivative
-    }
-    else // limited density
-    {
-        double a = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradUC * gradUC);
-        double beta = mInf * mInf * gradUC / a;
-        return -sqrt(a * a * a * a * a) / ((1 + beta * (gradU / gradUC - 1)) * (1 + beta * (gradU / gradUC - 1))) * beta / gradUC / gradU; // has to be multiplied by grad u to recover true derivative
-    }
+/**
+ * @brief Evaluate the isentropic density
+ *
+ * rho = (1 + (gamma-1)/2*Mi^2*(1-u^2))^(1/(gamma-1))
+ */
+double F0ElRho::eval(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+    double u = e.computeGradient(phi, k).norm(); // norm of velocity
+    double a2 = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u * u);
+    return sqrt(a2 * a2 * a2 * a2 * a2); // particularized to gamma=1.4
+}
+/**
+ * @brief Evaluate the isentropic density gradient
+ *
+ * drho = -Mi^2 * rho^(2-gamma) * u
+ * Computed value has to be multiplied by u to recover actual gradient
+ */
+double F0ElRho::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+    double u = e.computeGradient(phi, k).norm(); // norm of velocity
+    double a2 = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u * u);
+    return -mInf * mInf * sqrt(a2 * a2 * a2); // particularized to gamma=1.4
 }
 
 /**
- * @brief Evaluate the linear density (constant over element)
+ * @brief Evaluate the (clamped) isentropic density
+ *
+ * rho = a^(2/(gamma-1))
  */
-double F0ElRhoL::eval(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElRhoClamp::eval(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+    double a2 = F0ElSpeedSound::eval(e, phi, k);
+    return sqrt(a2 * a2 * a2 * a2 * a2); // particularized to gamma=1.4
+}
+/**
+ * @brief Evaluate the (clamped) isentropic density gradient
+ *
+ * drho = 1/(gamma-1) * a^(2*(2-gamma)) * da^2
+ * Computed value has to be multiplied by u to recover actual gradient
+ */
+double F0ElRhoClamp::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+    double a2 = F0ElSpeedSound::eval(e, phi, k);
+    double da2 = F0ElSpeedSound::evalGrad(e, phi, k);
+    return 1 / (gamma - 1) * sqrt(a2 * a2 * a2) * da2; // particularized to gamma=1.4
+}
+
+/**
+ * @brief Evaluate the linear density
+ */
+double F0ElRhoLin::eval(Element const &e, std::vector<double> const &phi, size_t k) const
 {
     return 1.;
 }
 /**
- * @brief Evaluate the linear density derivative (constant over element)
+ * @brief Evaluate the linear density gradient
  */
-double F0ElRhoL::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElRhoLin::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
 {
     return 0.;
 }
 
 /**
- * @brief Evaluate the nonlinear Mach number (constant over element)
+ * @brief Evaluate the isentropic Mach number
+ *
+ * M = u / sqrt(1/mi^2 + (gamma-1)/2*(1-u^2))
+ */
+double F0ElMach::eval(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+    double u2 = e.computeGradient(phi, k).squaredNorm();          // velocity squared
+    double a2 = 1 / (mInf * mInf) + 0.5 * (gamma - 1) * (1 - u2); // speed of sound squared
+    return sqrt(u2) / sqrt(a2);
+}
+/**
+ * @brief Evaluate the isentropic Mach number gradient
+ *
+ * dM = (1/(u*a) + (gamma-1)/2*u/a^3) * u
+ * Computed value has to be multiplied by u to recover actual gradient
  */
-double F0ElMach::eval(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElMach::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
 {
-    Eigen::VectorXd gradU = e.computeGradient(u, k);                  // velocity
-    double gradU2 = gradU.squaredNorm();                              // velocity squared
-    double a2 = 1 / (mInf * mInf) + 0.5 * (gamma - 1) * (1 - gradU2); // speed of sound squared
+    double u2 = e.computeGradient(phi, k).squaredNorm();          // velocity squared
+    double a2 = 1 / (mInf * mInf) + 0.5 * (gamma - 1) * (1 - u2); // speed of sound squared
+    return 1 / sqrt(u2 * a2) + 0.5 * (gamma - 1) * sqrt(u2) / sqrt(a2 * a2 * a2);
+}
 
-    return sqrt(gradU2) / sqrt(a2);
+/**
+ * @brief Evaluate the (clamped) isentropic Mach number
+ *
+ * M = u / sqrt(a^2/Mi^2)
+ */
+double F0ElMachClamp::eval(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+    double u = e.computeGradient(phi, k).norm(); // norm of velocity
+    double a2 = F0ElSpeedSound::eval(e, phi, k);
+    return mInf * u / sqrt(a2);
 }
 /**
- * @brief Evaluate the nonlinear Mach number derivative (constant over element)
+ * @brief Evaluate the (clamped) isentropic Mach number gradient
+ *
+ * dM = 1 / sqrt(a^2/Mi^2) + u / sqrt(a^6/Mi^2) * da^2
+ * Computed value has to be multiplied by u to recover actual gradient
  */
-double F0ElMach::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElMachClamp::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
 {
-    Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
-    double gradU2 = gradU.squaredNorm();             // velocity squared
-    double a2 = 1 / (mInf * mInf) + 0.5 * (gamma - 1) * (1 - gradU2);
-
-    return 1 / sqrt(gradU2 * a2) + 0.5 * (gamma - 1) * sqrt(gradU2) / sqrt(a2 * a2 * a2); // has to be multiplied by grad u to recover true derivative
+    double u = e.computeGradient(phi, k).norm(); // norm of velocity
+    double a2 = F0ElSpeedSound::eval(e, phi, k);
+    double da2 = F0ElSpeedSound::evalGrad(e, phi, k);
+    return mInf * (1 / sqrt(a2) / u - 0.5 * u / sqrt(a2 * a2 * a2) * da2);
 }
 
 /**
- * @brief Evaluate the linear Mach number (constant over element)
+ * @brief Evaluate the linear Mach number
  */
-double F0ElMachL::eval(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElMachLin::eval(Element const &e, std::vector<double> const &phi, size_t k) const
 {
     return 0.;
 }
 /**
- * @brief Evaluate the linear Mach number derivative (constant over element)
+ * @brief Evaluate the linear Mach number gradient
  */
-double F0ElMachL::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElMachLin::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
 {
     return 0.;
 }
 
 /**
- * @brief Evaluate the nonlinear pressure coefficient (constant over element)
+ * @brief Evaluate the isentropic pressure coefficient
+ *
+ * Cp = 2/(gamma*Mi^2) * ((1 + (gamma-1)/2*Mi^2*(1-u^2))^(gamma/(gamma-1)) - 1)
  */
-double F0ElCp::eval(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElCp::eval(Element const &e, std::vector<double> const &phi, size_t k) const
 {
-    Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
-    double gradU2 = gradU.squaredNorm();             // velocity squared
-    //particularized: 2 / (gamma * mInf * mInf) * (pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), gamma / (gamma - 1)) - 1);
-    double a = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2);
-
-    if (a <= 0)
-        return -2 / (gamma * mInf * mInf); // limit the pressure coefficient when vacuum conditions are reached
-    else
-        return 2 / (gamma * mInf * mInf) * (sqrt(a * a * a * a * a * a * a) - 1);
+    double u2 = e.computeGradient(phi, k).squaredNorm(); // velocity squared
+    double a2 = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u2);
+    return 2 / (gamma * mInf * mInf) * (sqrt(a2 * a2 * a2 * a2 * a2 * a2 * a2) - 1); // particularized to gamma=1.4
 }
 /**
- * @brief Evaluate the nonlinear pressure coefficient derivative (constant over element)
+ * @brief Evaluate the isentropic pressure coefficient gradient
+ *
+ * dCp = -2*(1 + (gamma-1)/2*Mi^2*(1-u^2))^(1/(gamma-1)) * u
+ * Computed value has to multiplied by u to recover actual gradient
  */
-double F0ElCp::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElCp::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
 {
-    Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
-    double gradU2 = gradU.squaredNorm();             // velocity squared
-    //particularized: -2 * pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), 1 / (gamma - 1));
-    double a = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2);
+    double u2 = e.computeGradient(phi, k).squaredNorm(); // velocity squared
+    double a2 = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u2);
+    return -2 * sqrt(a2 * a2 * a2 * a2 * a2); // particularized to gamma=1.4
+}
 
-    if (a <= 0)
-        return 0; // limit the gradient when vacuum conditions are reached
-    else
-        return -2 * sqrt(a * a * a * a * a); // has to be multiplied by grad u to recover true derivative
+/**
+ * @brief Evaluate the (clamped) isentropic pressure coefficient
+ *
+ * Cp = 2/(gamma*Mi^2) * (a^(2*gamma/(gamma-1)) - 1)
+ */
+double F0ElCpClamp::eval(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+    double a2 = F0ElSpeedSound::eval(e, phi, k);
+    return 2 / (gamma * mInf * mInf) * (sqrt(a2 * a2 * a2 * a2 * a2 * a2 * a2) - 1); // particularized to gamma=1.4
+}
+/**
+ * @brief Evaluate the (clamped) isentropic pressure coefficient gradient
+ *
+ * dCp = 2/((gamma-1)*Mi^2) * a^(2/(gamma-1)) * da^2
+ * Computed value has to multiplied by u to recover actual gradient
+ */
+double F0ElCpClamp::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
+{
+    double a2 = F0ElSpeedSound::eval(e, phi, k);
+    double da2 = F0ElSpeedSound::evalGrad(e, phi, k);
+    return 2 / ((gamma - 1) * mInf * mInf) * sqrt(a2 * a2 * a2 * a2 * a2) * da2; // particularized to gamma=1.4
 }
 
 /**
- * @brief Evaluate the linear pressure coefficient (constant over element)
+ * @brief Evaluate the linear pressure coefficient
+ *
+ * Cp = 1-u^2
  */
-double F0ElCpL::eval(Element const &e, std::vector<double> const &u, size_t k) const
+double F0ElCpLin::eval(Element const &e, std::vector<double> const &phi, size_t k) const
 {
-    Eigen::VectorXd gradU = e.computeGradient(u, k); // velocity
-    double gradU2 = gradU.squaredNorm();             // velocity squared
-    return 1 - gradU2;
+    return 1 - e.computeGradient(phi, k).squaredNorm();
 }
-double F0ElCpL::evalGrad(Element const &e, std::vector<double> const &u, size_t k) const
+/**
+ * @brief Evaluate the linear pressure coefficient gradient
+ *
+ * dCp = -2*u
+ * Computed value has to multiplied by u to recover actual gradient
+ */
+double F0ElCpLin::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
 {
-    return -2.; // has to be multiplied by grad u to recover true derivative
+    return -2.;
 }
\ No newline at end of file

dart/src/wF0El.h

@@ -34,9 +34,9 @@ public:
     F0El() : Observer() {}
     virtual ~F0El() {}
 
-    virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+    virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
                         size_t k) const = 0;
-    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
                             size_t k) const = 0;
 };
 
@@ -51,55 +51,96 @@ public:
     F0ElC(double _v) : F0El(), v(_v) {}
     virtual void update(double _v) override;
 
-    virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+    virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
                         size_t k) const override;
-    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
                             size_t k) const override;
 };
 
 /**
- * @brief Nonlinear density
+ * @brief Speed of sound squared
  *
- * The density is limited (Padé approximation) for large value of velocity and
- * particularized for diatomic gas. The input _mC2 is the square of the critical
- * Mach number used to compute the critical velocity.
+ * This function computes the square of the speed of sound, and is then used
+ * to compute the actual density, Mach number and pressure coefficient. This
+ * function is clamped for large velocities using a Padé approximation and
+ * particularized to diatomic gas. The function is normalized by its
+ * freestream value.
+ */
+class DART_API F0ElSpeedSound : public F0El
+{
+protected:
+    double gamma; ///< heat capacity ratio (diatomic gas only)
+    double mInf;  ///< freestream Mach number
+private:
+    double uC;   ///< critical velocity
+    double aC;   ///< critical speed of sound squared
+    double beta; ///< ratio for Padé
+
+public:
+    F0ElSpeedSound(double _mInf, double _mC2 = 5);
+
+    virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
+                        size_t k) const override;
+    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
+                            size_t k) const override;
+};
+
+/**
+ * @brief Isentropic density
+ *
+ * The density is particularized to diatomic gas. The density is normalized by
+ * its freestream value.
  */
 class DART_API F0ElRho : public F0El
 {
-    double gamma;  ///< heat capacity ratio (diatomic gas only)
-    double mInf;   ///< freestream Mach number
-    double gradUC; ///< critical velocity
+    double gamma; ///< heat capacity ratio (diatomic gas only)
+    double mInf;  ///< freestream Mach number
 
 public:
-    F0ElRho(double _mInf, double _mC2 = 5) : F0El(), gamma(1.4), mInf(_mInf)
-    {
-        gradUC = sqrt(_mC2 * (1 / (mInf * mInf) + (gamma - 1) / 2) / (1 + (gamma - 1) / 2 * _mC2));
-    }
+    F0ElRho(double _mInf) : F0El(), gamma(1.4), mInf(_mInf) {}
+
+    virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
+                        size_t k) const override;
+    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
+                            size_t k) const override;
+};
+
+/**
+ * @brief Isentropic clamped density
+ *
+ * The density is clamped (Padé approximation) for large values of velocity and
+ * particularized to diatomic gas. The density is normalized by its
+ * freestream value.
+ */
+class DART_API F0ElRhoClamp : public F0ElSpeedSound
+{
+public:
+    F0ElRhoClamp(double _mInf, double _mC2 = 5) : F0ElSpeedSound(_mInf, _mC2) {}
 
-    virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+    virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
                         size_t k) const override;
-    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
                             size_t k) const override;
 };
 
 /**
  * @brief Linear density (constant)
  */
-class DART_API F0ElRhoL : public F0El
+class DART_API F0ElRhoLin : public F0El
 {
 public:
-    F0ElRhoL() : F0El() {}
+    F0ElRhoLin() : F0El() {}
 
-    virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+    virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
                         size_t k) const override;
-    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
                             size_t k) const override;
 };
 
 /**
- * @brief Nonlinear Mach number
+ * @brief Isentropic Mach number
  *
- * Particularized for diatomic gas
+ * The Mach number is particularized to diatomic gas.
  */
 class DART_API F0ElMach : public F0El
 {
@@ -109,30 +150,47 @@ class DART_API F0ElMach : public F0El
 public:
     F0ElMach(double _mInf) : F0El(), gamma(1.4), mInf(_mInf) {}
 
-    virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+    virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
+                        size_t k) const override;
+    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
+                            size_t k) const override;
+};
+
+/**
+ * @brief Isentropic clamped Mach number
+ *
+ * The Mach number is clamped (Padé approximation) for large values of velocity
+ * and particularized to diatomic gas.
+ */
+class DART_API F0ElMachClamp : public F0ElSpeedSound
+{
+public:
+    F0ElMachClamp(double _mInf, double _mC2 = 5) : F0ElSpeedSound(_mInf, _mC2) {}
+
+    virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
                         size_t k) const override;
-    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
                             size_t k) const override;
 };
 
 /**
  * @brief Linear Mach number (constant)
  */
-class DART_API F0ElMachL : public F0El
+class DART_API F0ElMachLin : public F0El
 {
 public:
-    F0ElMachL() : F0El() {}
+    F0ElMachLin() : F0El() {}
 
-    virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+    virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
                         size_t k) const override;
-    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
                             size_t k) const override;
 };
 
 /**
- * @brief Nonlinear pressure coefficient
+ * @brief Isentropic pressure coefficient
  *
- * Particularized for diatomic gas
+ * The pressure coefficient is particularized to diatomic gas.
  */
 class DART_API F0ElCp : public F0El
 {
@@ -142,27 +200,43 @@ class DART_API F0ElCp : public F0El
 public:
     F0ElCp(double _mInf) : F0El(), gamma(1.4), mInf(_mInf) {}
 
-    virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+    virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
                         size_t k) const override;
-    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
                             size_t k) const override;
 };
 
 /**
- * @brief Linear pressure coefficient
+ * @brief Isentropic clamped pressure coefficient
+ *
+ * The pressure coefficient is clamped (Padé approximation) for large values of
+ * velocity and particularized to diatomic gas.
  */
-class DART_API F0ElCpL : public F0El
+class DART_API F0ElCpClamp : public F0ElSpeedSound
 {
+public:
+    F0ElCpClamp(double _mInf, double _mC2 = 5) : F0ElSpeedSound(_mInf, _mC2) {}
 
+    virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
+                        size_t k) const override;
+    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
+                            size_t k) const override;
+};
+
+/**
+ * @brief Linear pressure coefficient
+ */
+class DART_API F0ElCpLin : public F0El
+{
 public:
-    F0ElCpL() : F0El() {}
+    F0ElCpLin() : F0El() {}
 
-    virtual double eval(tbox::Element const &e, std::vector<double> const &u,
+    virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
                         size_t k) const override;
-    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &u,
+    virtual double evalGrad(tbox::Element const &e, std::vector<double> const &phi,
                             size_t k) const override;
 };
 
 } // namespace dart
 
-#endif //WF0EL_H
+#endif // WF0EL_H

dart/src/wF0Ps.h

@@ -54,4 +54,4 @@ public:
 
 } // namespace dart
 
-#endif //WF0PS_H
+#endif // WF0PS_H

dart/src/wF1Ct.h

@@ -50,7 +50,7 @@ class DART_API F1CtDrag : public F1Ct
 
 public:
     F1CtDrag(int _nDim, double _sRef, double alpha, double _beta = 0);
-    void update(double alpha);
+    virtual void update(double alpha) override;
 
     virtual Eigen::Vector3d eval() const override;
     virtual Eigen::Vector3d evalGrad() const override;
@@ -69,7 +69,7 @@ class DART_API F1CtSide : public F1Ct
 
 public:
     F1CtSide(int _nDim, double _sRef, double alpha, double _beta = 0);
-    void update(double alpha);
+    virtual void update(double alpha) override;
 
     virtual Eigen::Vector3d eval() const override;
     virtual Eigen::Vector3d evalGrad() const override;
@@ -88,7 +88,7 @@ class DART_API F1CtLift : public F1Ct
 
 public:
     F1CtLift(int _nDim, double _sRef, double alpha, double _beta = 0);
-    void update(double alpha);
+    virtual void update(double alpha) override;
 
     virtual Eigen::Vector3d eval() const override;
     virtual Eigen::Vector3d evalGrad() const override;
@@ -96,4 +96,4 @@ public:
 
 } // namespace dart
 
-#endif //WF1CT_H
+#endif // WF1CT_H

dart/src/wF1El.h

@@ -62,4 +62,4 @@ public:
 
 } // namespace dart
 
-#endif //WF1EL_H
+#endif // WF1EL_H

dart/src/wFace.h

@@ -85,4 +85,4 @@ public:
 
 } // namespace dart
 
-#endif //WFACE_H
+#endif // WFACE_H

dart/src/wFluid.cpp

@@ -54,15 +54,15 @@ void Fluid::initFun(double mInf, int dim, double alpha, double beta)
 {
     if (mInf == 0)
     {
-        rho = new F0ElRhoL();
-        mach = new F0ElMachL();
-        cP = new F0ElCpL();
+        rho = new F0ElRhoLin();
+        mach = new F0ElMachLin();
+        cP = new F0ElCpLin();
     }
     else
     {
-        rho = new F0ElRho(mInf);
-        mach = new F0ElMach(mInf);
-        cP = new F0ElCp(mInf);
+        rho = new F0ElRhoClamp(mInf, 3.0);
+        mach = new F0ElMachClamp(mInf, 3.0);
+        cP = new F0ElCpClamp(mInf, 3.0);
     }
     phiInf = new F0PsPhiInf(dim, alpha, beta);
 }

dart/src/wFluid.h

@@ -59,4 +59,4 @@ public:
 
 } // namespace dart
 
-#endif //WFLUID_H
+#endif // WFLUID_H

dart/src/wFpeLSFunction.h

@@ -49,4 +49,4 @@ public:
 
 } // namespace dart
 
-#endif //WFPELSFUNCTION_H
+#endif // WFPELSFUNCTION_H

dart/src/wFreestream.h

@@ -45,4 +45,4 @@ public:
 
 } // namespace dart
 
-#endif //WFREESTREAM_H
+#endif // WFREESTREAM_H

dart/src/wFreestreamResidual.h

@@ -38,4 +38,4 @@ public:
 };
 
 } // namespace dart
-#endif //WFREESTREAMRESIDUAL_H
+#endif // WFREESTREAMRESIDUAL_H

dart/src/wKutta.cpp

@@ -82,10 +82,10 @@ void Kutta::connectNodes()
     // Create the node map
     for (auto nUp : nodesUp)
     {
-        //std::cout << "MASTER: " << nUp->pos(0) << ',' << nUp->pos(1) << ',' << nUp->pos(2) << std::endl;
+        // std::cout << "MASTER: " << nUp->pos(0) << ',' << nUp->pos(1) << ',' << nUp->pos(2) << std::endl;
         for (auto nLw : nodesLw)
         {
-            //std::cout << "    SLAVE: " << nLw->pos(0) << ',' << nLw->pos(1) << ',' << nLw->pos(2) << std::endl;
+            // std::cout << "    SLAVE: " << nLw->pos(0) << ',' << nLw->pos(1) << ',' << nLw->pos(2) << std::endl;
             if ((nUp->pos - nLw->pos).norm() <= 1e-14)
             {
                 nodMap.push_back(std::pair<Node *, Node *>(nUp, nLw));

dart/src/wKutta.h

@@ -56,4 +56,4 @@ private:
 
 } // namespace dart
 
-#endif //WKUTTA_H
+#endif // WKUTTA_H