Handled degenrate case in Kutta

hspm/src/geometry.cpp

@@ -30,7 +30,7 @@ void HSPM::generateNaca4DigitCoordinates(double camber, double camberPos, double
 		} else {
 			x_distr(i) = chord/2 * (cos(angle)+1);
 		}*/
-		x_distr(i) = chord * pow((1 + cos((angle + M_PI) / 2)),2);
+		x_distr(i) = chord * pow((1 + cos((angle + M_PI) / 2)), 2);
 	}
 
 	for (size_t i=0; i<=N; i++) {

hspm/src/solver.cpp

@@ -34,7 +34,7 @@ void HSPM::solve()
     // q = s1*tau + s2
     // We need Kutta for tau
     tau = this->solveOffBodyKutta();
-    
+
     q = s1*tau + s2;
 
     this->computeInviscidVelocity(); 
@@ -84,6 +84,9 @@ double HSPM::solveOffBodyKutta() {
     double _tau_terms = 2*(_a*_b + _c*_d - _e*_f - _g*_h);
     double _const_terms = pow(_b,2) + pow(_d,2) - pow(_f,2) - pow(_h,2);
 
+    if (_tau2_terms < 1e-10) 
+        return -_const_terms / _tau_terms; // If the quadratic equation is degenerate, we can solve it directly
+
     // Solve the quadratic equation
     double tau_ = (-_tau_terms + sqrt(pow(_tau_terms,2) - 4*_tau2_terms*_const_terms)) / (2*_tau2_terms);
 

python/main.py

 from utils import *
 
-N = 1000
-
 config = {
 	'chord': 1,
-	'it_solver_tolerance': 1e-6,
 	'aoa': 2,
 	'naca': "0012",
-	'N': N,
+	'N': 100,
 	'filename': "airfoil.dat"
 }
 
 if __name__ == "__main__":
 	solver = initHSPM(config)
 
-	"""for i in range(N):
-					solver.imposeBlowingVelocity(i, 1)
-					solver.setdStar(i, 0.1)"""
-
 	solver.solve()
 
 	cl = solver.cl