diff --git a/hspm/api/core.py b/hspm/api/core.py
index 53a19be82ed60bebdd056bcde83b406c1357411a..fd21cc302f6579ab00f3e817118e4ca84e027c12 100644
--- a/hspm/api/core.py
+++ b/hspm/api/core.py
@@ -1,6 +1,30 @@
+mesh_types = {
+ "CosN": 0,
+ "ConstantN": 1,
+ "Cos2N": 2,
+ "LinearS": 3,
+ "ExponentialS": 4
+}
+
def initHSPM(cfg):
"""
Initializes the HSPM object with the given configuration
+
+ Parameters
+ ----------
+ cfg : dict
+ Configuration dictionary
+ Contains:
+ - chord (float): Chord length
+ - aoa (float): Angle of attack
+
+ * MESH PARAMETERS
+ - meshType (str): Type of mesh (CosN, ConstantN, Cos2N, LinearS, ExponentialS)
+ - sLe (float): Leading edge mesh size, only for S-type meshes
+ - sTe (float): Trailing edge mesh size, only for S-type meshes
+ - naca (str): NACA airfoil code
+ or
+ - filename (str): Airfoil file names
"""
import hspmw
@@ -8,10 +32,28 @@ def initHSPM(cfg):
hspm = hspmw.HSPM()
- hspm.chord = cfg['chord'] # Not optimal, but it will do for now
- hspm.V_inf = 1 # Shouldn't have any impact on the solution
+ hspm.chord = cfg['chord']
+ hspm.V_inf = 1
+ if not 'aoa' in cfg:
+ print("HSPM Warning: aoa not defined. Using default value of 0")
+ cfg['aoa'] = 0
hspm.AoA = np.deg2rad(cfg['aoa'])
+ if not "meshType" in cfg:
+ print("HSPM Warning: meshType not defined. Using default cosine mesh (CosN)")
+ cfg['meshType'] = "CosN"
+ hspm.meshType = mesh_types[cfg['meshType']]
+ if cfg['meshType'] == "LinearS" or cfg['meshType'] == "ExponentialS":
+ if not "sLe" in cfg:
+ print("HSPM Warning: Leading edge mesh size (sLe) not defined. Using default value of 0.01")
+ cfg['sLe'] = 0.01
+ if not "sTe" in cfg:
+ print("HSPM Warning: Trailing edge mesh size (sTe) not defined. Using default value of 0.01")
+ cfg['sTe'] = 0.01
+ hspm.sLe = cfg['sLe']
+ hspm.sTe = cfg['sTe']
+ cfg["N"] = 1 # Placeholder value for the generateNaca4DigitCoordinates function
+
if cfg['naca']:
hspm.generateNaca4DigitCoordinates(int(cfg['naca'][0]), int(cfg['naca'][1]), int(cfg['naca'][2:4]), int(cfg['N']))
else:
diff --git a/hspm/src/geometry.cpp b/hspm/src/geometry.cpp
index efd5b934bfa0aea68fc85d40bd73c1350be1724c..6c1321e12a42af964317dde7aa49e6ba7de8b3ba 100644
--- a/hspm/src/geometry.cpp
+++ b/hspm/src/geometry.cpp
@@ -1,4 +1,5 @@
#include "geometry.h"
+#include <iostream>
void HSPM::generateNaca4DigitCoordinates(double camber, double camberPos, double thickness, size_t _N)
{
@@ -14,23 +15,106 @@ void HSPM::generateNaca4DigitCoordinates(double camber, double camberPos, double
Thickness of the airfoil [%]. Third/fourth digits
*/
+ if (meshType == LinearS) {
+ // Figure out how many elements are required
+ double N_double = 4.0/(sLe + sTe);
+ // Round to the nearest even number
+ _N = (size_t) ((int) N_double/2.0) * 2;
+ std::cout << "HSPM: The number of elements has been set to " << _N << " to match the prescribed sizes at the LE and TE." << std::endl;
+ } else if (meshType == ExponentialS)
+ {
+ // Figure out how many elements are required
+ if (sTe == sLe) { // In this case, we will use constant size elements (might be slightly off the prescribed size, but minimal variation is expected)
+ double N_double = 2.0 / sLe;
+
+ // Round to the nearest even number
+ _N = (size_t) ((int) N_double/2.0) * 2;
+ std::cout << "HSPM: The number of elements has been set to " << _N << " to match the prescribed sizes at the LE and TE." << std::endl;
+
+ sLe = 2.0 / _N;
+ sTe = 2.0 / _N;
+
+ prog = 1.0;
+ }
+ else {
+ // Get the progression
+ double _r = sTe - sLe + 1;
+ double N_double = 2.0 * (log(sTe/sLe) / log(_r) + 1.0);
+ // Round to the nearest even number
+ _N = (size_t) ((int) N_double/2.0) * 2;
+ std::cout << "HSPM: The number of elements has been set to " << _N << " to match the prescribed sizes at the LE and TE." << std::endl;
+
+ // Compute the value of the progression (needs to solve a nonlinear equation w/ Newton's method)
+ // 1-r = sLe * (1-r^(N/2))
+ auto equation = [this](double r, double _N) {
+ return 1-r - sLe * (1-pow(r, _N/2.0));
+ };
+
+ auto equationDerivative = [this](double r, double _N) {
+ return -1 - sLe * (-pow(r, _N/2.0-1) * _N/2.0);
+ };
+
+ double r0 = _r;
+ for (int i = 0; i < 1000; ++i) {
+ double h = equation(r0, _N) / equationDerivative(r0, _N);
+ if (std::abs(h) < 1e-8) {
+ break;
+ }
+ r0 = r0 - h;
+ }
+ prog = r0;
+ }
+ }
+
N = _N;
x = Eigen::VectorXd(N+1);
y = Eigen::VectorXd(N+1);
Eigen::VectorXd angleDistribution = Eigen::VectorXd::LinSpaced(N+1, 0, 2*M_PI);
-
+
Eigen::VectorXd x_distr = Eigen::VectorXd::Zero(N+1);
for (size_t i=0; i<=N; i++) {
double angle = angleDistribution(i);
- /*if (angle < M_PI/2) {
- x_distr(i) = chord/2 * (1 - angle * 2 / M_PI + 1);
- } else if (angle > 3*M_PI/2) {
- x_distr(i) = chord/2 * ((angle - 3*M_PI/2) * 2 / M_PI + 1);
- } else {
- x_distr(i) = chord/2 * (cos(angle)+1);
- }*/
- x_distr(i) = chord * pow((1 + cos((angle + M_PI) / 2)), 2);
+
+ if (meshType == CosN) // Cosine mesh with N elements
+ x_distr(i) = chord/2 * (1 + cos(angle));
+
+ else if (meshType == ConstantN) // Constant size mesh with N elements
+ {
+ if (angle < M_PI)
+ x_distr(i) = chord * (-M_1_PI * angle + 1);
+ else x_distr(i) = chord * (M_1_PI * angle - 1);
+ }
+
+ else if (meshType == Cos2N) // Cosine square mesh with N elements
+ x_distr(i) = chord * pow((1 + cos((angle + M_PI) / 2)), 2);
+
+ else if (meshType == LinearS) // Linear size mesh, with sizes at LE and TE (approx) prescribed
+ { // Maths au milieu de mon cahier
+ // Compute the size increment (could be saved to memory, but mesh generation is fast enough, so it's fine here)
+ double delta_s = (1.0 - (N/2.0)*sLe) / ((N/2.0*(N/2.0-1.0))/2.0);
+ if (i <= N/2.0) { // On the top side, from LE to TE
+ size_t j = N/2.0-i;
+
+ x_distr(j) = sLe * i + delta_s * (i-1) * i / 2.0;
+
+ } else {
+ x_distr(i) = x_distr(i-1) + sLe + delta_s * (i-(N/2.0+1));
+ }
+ }
+
+ else if (meshType == ExponentialS) // Exponential size mesh, with sizes @ LE and TE (approx) prescribed
+ {
+ if (i <= N/2.0) { // Top side
+ size_t j = N/2.0-i;
+
+ if (prog == 1) x_distr(j) = sLe * i;
+ else x_distr(j) = sLe * (pow(prog, i) - 1) / (prog - 1);
+
+ } else { // Bottom side
+ x_distr(i) = x_distr(i-1) + sLe * pow(prog, i-(N/2.0+1));
+ }
+ }
}
for (size_t i=0; i<=N; i++) {
@@ -97,6 +181,8 @@ void HSPM::loadCoordinates(char* fileName) {
N = coordinates.cols() - 1;
x = coordinates.row(0);
y = coordinates.row(1);
+
+ // TODO: will need to close some airfoils if the first and last points are not the same !!!
}
void HSPM::saveCoordinates(char* fileName) {
diff --git a/hspm/src/hspm.h b/hspm/src/hspm.h
index 640c1ba1aa8cab65c79c45c2ec8a4ce05bb08299..0e37ce4df2414947e8535b2de8ecf31e18cd5c8e 100644
--- a/hspm/src/hspm.h
+++ b/hspm/src/hspm.h
@@ -71,6 +71,19 @@ public:
Eigen::VectorXd y_m_wake;
Eigen::VectorXd U_wake;
Eigen::VectorXd V_wake;
+
+ // Mesh stuff
+ enum MeshType {
+ CosN = 0,
+ ConstantN = 1,
+ Cos2N = 2,
+ LinearS = 3,
+ ExponentialS = 4
+ };
+ MeshType meshType;
+ double sLe, sTe;
+ double prog; // For exponential mesh
+
private:
Eigen::PartialPivLU<Eigen::MatrixXd> solver;
};