From 967bfb1edc33eb60489f5bb2011cb42c20b07764 Mon Sep 17 00:00:00 2001
From: Paul <paul.dechamps@student.uliege.be>
Date: Sun, 10 Jul 2022 17:26:17 +0200
Subject: [PATCH] Using 2 meshes for 3D
---
sdpm/src/pPanels.cpp | 23 +-
sdpm/src/pPanels.h | 21 +-
sdpm/src/pSDSolver.cpp | 74 ++-
sdpm/src/pSDSolver.h | 3 +-
sdpm/tests/panels3d.py | 6 +-
vii/models/n0012_3_visc.geo | 397 +++++++++++++
vii/pyVII/vInterpolator.py | 1053 +++++++++++++++++++++++------------
vii/pyVII/vUtils.py | 24 +
vii/src/wBLRegion.cpp | 42 +-
vii/src/wClosures.cpp | 119 +++-
vii/src/wClosures.h | 5 +-
vii/src/wTimeSolver.cpp | 6 +-
vii/src/wViscFluxes.cpp | 6 +-
vii/src/wViscSolver.cpp | 11 +-
vii/tests/bli2.py | 2 +-
vii/tests/bli3.py | 237 ++++----
16 files changed, 1436 insertions(+), 593 deletions(-)
create mode 100644 vii/models/n0012_3_visc.geo
diff --git a/sdpm/src/pPanels.cpp b/sdpm/src/pPanels.cpp
index f930572..c565b9c 100644
--- a/sdpm/src/pPanels.cpp
+++ b/sdpm/src/pPanels.cpp
@@ -9,12 +9,15 @@ using namespace sdpm;
Panels::Panels(size_t _nPanChord)
{
nPanChord = _nPanChord;
- nPanSpan = 2;
+ nPanSpan = 20;
beta = 1;
}
void Panels::CreateWing()
{
+
+ b = AR*(c0+c0*_lambda)/2;
+
Eigen::RowVectorXd xp(2*nPanChord+1); // Non-linear chordwise distribution
Eigen::RowVectorXd xc(2*nPanChord+1); // Control point chordwise distribution
@@ -46,7 +49,7 @@ void Panels::CreateWing()
/* Spanwise chord variation */
Eigen::RowVectorXd cvar;
- cvar = (lambda-1)*c0*yp.cwiseAbs();
+ cvar = (_lambda-1)*c0*yp.cwiseAbs();
cvar += c0 * Eigen::MatrixXd::Ones(1, cvar.size());
Eigen::RowVectorXd epsilonvar;
epsilonvar=epsilon*yp.cwiseAbs();
@@ -63,26 +66,24 @@ void Panels::CreateWing()
Eigen::MatrixXd zp = Eigen::MatrixXd::Zero(1,2*nPanChord+1);
- double t = .12; // Airfoil maximum thickness
-
/* Upper side */
for (auto i=0; i<zp.size()/2; ++i)
{
- zp(i) = -(t/0.2)*(0.2969*std::sqrt(xp(i)) - .126 * xp(i) - .3516* xp(i)*xp(i) + .2843*xp(i)*xp(i)*xp(i) - .1036*xp(i)*xp(i)*xp(i)*xp(i));
+ zp(i) = -(thick/0.2)*(0.2969*std::sqrt(xp(i)) - .126 * xp(i) - .3516* xp(i)*xp(i) + .2843*xp(i)*xp(i)*xp(i) - .1036*xp(i)*xp(i)*xp(i)*xp(i));
}
/* Lower side */
for (auto i=zp.size()/2; i<zp.size(); ++i)
{
- zp(i) = (t/0.2)*(0.2969*std::sqrt(xp(i)) - .126 * xp(i) - .3516* xp(i)*xp(i) + .2843*xp(i)*xp(i)*xp(i) - .1036*xp(i)*xp(i)*xp(i)*xp(i));
+ zp(i) = (thick/0.2)*(0.2969*std::sqrt(xp(i)) - .126 * xp(i) - .3516* xp(i)*xp(i) + .2843*xp(i)*xp(i)*xp(i) - .1036*xp(i)*xp(i)*xp(i)*xp(i));
}
- // iko=find(xp <= maxpos);
- // zp(iko)=zp(iko)+maxcamb/maxpos^2*(2*maxpos*xp(iko)-xp(iko).^2);
- // iko=find(xp > maxpos);
- // zp(iko)=zp(iko)+maxcamb/(1-maxpos)^2*((1-2*maxpos)+2*maxpos*xp(iko)-xp(iko).^2);
+ /* iko = xp <= maxpos;
+ zp(iko)=zp(iko)+maxcamb/(maxpos*maxpos)*(2*maxpos*xp(iko)-xp(iko).array()*xp(iko).array());
+ iko=xp > maxpos;
+ zp(iko)=zp(iko)+maxcamb/((1-maxpos)*(1-maxpos))*((1-2*maxpos)+2*maxpos*xp(iko)-xp(iko).array()*xp(iko).array()); */
/* Assemble bound panel corner point matrices */
/* Eigen::Array * Eigen::Array is a mult element wise while Eigen::Matrix * Eigen::Matrix is a matrix product -> cwiseProduct */
@@ -90,7 +91,7 @@ void Panels::CreateWing()
Yp = yp.replicate(2*nPanChord+1,1) * b/2;
Zp = zp.transpose().replicate(1,nPanSpan+1).cwiseProduct(cvar.replicate(2*nPanChord+1,1));
- Surface = (c0+c0*lambda)/2 * b;
+ Surface = (c0+c0*_lambda)/2 * b;
diff --git a/sdpm/src/pPanels.h b/sdpm/src/pPanels.h
index 736330c..d613b65 100644
--- a/sdpm/src/pPanels.h
+++ b/sdpm/src/pPanels.h
@@ -3,6 +3,7 @@
#include "sdpm.h"
#include <vector>
+#include <cmath>
#include <Eigen/Dense>
namespace sdpm
@@ -14,13 +15,19 @@ class DART_API Panels
private:
public:
- double c0 = 1; /* Root chord (m) */
- double lambda = 1; /* Taper ratio (-) */
- double epsilon = 0; /* Twist angle (rad) */
- double LambdaC4 = 0; /* Sweep angle at 1/4 chord (rad) */
- double b = 1; /* Span (m) */
- double dihedral = 0; /* Dihedral angle (rad) */
- double ax = 0; /* Axis around which twist and taper are defined */
+ double c0 = 2; /* Root chord (m) */
+ double _lambda = 0.3; /* Taper ratio (-) */
+ double epsilon = 3*M_PI/180; /* Twist angle (rad) */
+ double LambdaC4 = 35*M_PI/180; /* Sweep angle at 1/4 chord (rad) */
+ double AR = 7;
+ double b; /* Span (m) */
+ double dihedral = 3*M_PI/180; /* Dihedral angle (rad) */
+ double ax = 1/2; /* Axis around which twist and taper are defined */
+
+ double maxcamb = 0.02;
+ double maxpos = 0.4;
+ double thick = .12; // Airfoil maximum thickness
+
double beta;
double Surface;
diff --git a/sdpm/src/pSDSolver.cpp b/sdpm/src/pSDSolver.cpp
index 721fc44..581108a 100644
--- a/sdpm/src/pSDSolver.cpp
+++ b/sdpm/src/pSDSolver.cpp
@@ -4,20 +4,28 @@
#include <vector>
#include <cmath>
+#include <iomanip>
+#include <deque>
+
using namespace sdpm;
-SDSolver::SDSolver(Panels *_pan, double _alpha, double _airspeed, unsigned long _maxIt)
+SDSolver::SDSolver(Panels *_pan, double _alpha, double _airspeed, double _tol, unsigned long _maxIt)
{
pan = _pan;
fs_alpha = _alpha;
fs_airspeed = _airspeed;
+ tol = _tol;
maxIt = _maxIt;
}
void SDSolver::Run()
{
- int R = 1; /* Wake type */
+ std::cout << std::fixed << std::setprecision(2)
+ << "\nSolving flow for alpha " << fs_alpha * 180 / 3.14159 << "deg AoA, "
+ << std::endl;
+
+ int R = 0; /* Wake type */
Eigen::Vector3d fs_velocity;
@@ -51,6 +59,10 @@ void SDSolver::Run()
pan->Yp0 = pan->Yp;
pan->Zp0 = pan->Zp;
+ Eigen::MatrixXd uc;
+ Eigen::MatrixXd vc;
+ Eigen::MatrixXd wc;
+
/* std::cout << " xp yp zp " << std::endl;
std::cout << pan->Xp << std::endl;
std::cout << pan->Yp << std::endl;
@@ -229,9 +241,17 @@ void SDSolver::Run()
unsigned long it = 1;
- while(it < maxIt){
+ std::cout << std::setw(8) << "Iters"
+ << std::setw(12) << "Cl"
+ << std::setw(12) << "Cd"
+ << std::setw(15) << "Error\n";
+ std::cout << std::fixed << std::setprecision(5);
+ std::cout << std::setw(8) << 0
+ << std::setw(12) << 0
+ << std::setw(12) << 1
+ << std::setw(15) << "-" << "\n";
- std::cout << it << std::endl;
+ while(it < maxIt){
/* Propagate wake */
Eigen::MatrixXd uwp = Eigen::MatrixXd::Zero(it-1, pan->nPanSpan + 1);
@@ -408,9 +428,9 @@ void SDSolver::Run()
/* Calculate velocities on control points in the x,y,z directions
On each control point we project the velocities in the chordwise,
spanwise and normal directions to the x,y,z axes */
- Eigen::MatrixXd uc(2*pan->nPanChord, pan->nPanSpan);
- Eigen::MatrixXd vc(2*pan->nPanChord, pan->nPanSpan);
- Eigen::MatrixXd wc(2*pan->nPanChord, pan->nPanSpan);
+ uc.setZero(2*pan->nPanChord, pan->nPanSpan);
+ vc.setZero(2*pan->nPanChord, pan->nPanSpan);
+ wc.setZero(2*pan->nPanChord, pan->nPanSpan);
Eigen::Matrix3d Rmat_sys;
Eigen::Vector3d RHS_sys;
@@ -443,45 +463,40 @@ void SDSolver::Run()
}
}
- uc = Eigen::MatrixXd::Constant(uc.rows(), uc.cols(), fs_velocity(0)) + uc;
- vc = Eigen::MatrixXd::Constant(vc.rows(), vc.cols(), fs_velocity(1)) + vc;
- wc = Eigen::MatrixXd::Constant(wc.rows(), wc.cols(), fs_velocity(2)) + wc;
+ uc += Eigen::MatrixXd::Constant(uc.rows(), uc.cols(), fs_velocity(0));
+ vc += Eigen::MatrixXd::Constant(vc.rows(), vc.cols(), fs_velocity(1));
+ wc += Eigen::MatrixXd::Constant(wc.rows(), wc.cols(), fs_velocity(2));
/* Calculate presssure coefficient */
Eigen::MatrixXd dPhidt;
if (it==1){
- Eigen::Map<Eigen::VectorXd> Phi1Flatten(Phi[1].data(),Phi[1].size());
- Eigen::Map<Eigen::VectorXd> Phi0Flatten(Phi0.data(),Phi0.size());
- dPhidt = (Phi1Flatten-Phi0Flatten)/dt;
+ dPhidt = (Phi[1]-Phi0)/dt;
}
else{
- Eigen::Map<Eigen::VectorXd> Phi1Flatten(Phi[it].data(),Phi[it].size());
- Eigen::Map<Eigen::VectorXd> Phi0Flatten(Phi[it-1].data(),Phi[it-1].size());
- dPhidt = (Phi1Flatten - Phi0Flatten)/dt;
+ dPhidt = (Phi[it] - Phi[it-1])/dt;
}
Eigen::MatrixXd cp_dum = (uc.array()*uc.array() + vc.array()*vc.array() + wc.array()*wc.array())*1/fs_airspeed/fs_airspeed - 2*dPhidt.array()*1/fs_airspeed/fs_airspeed;
cp = Eigen::MatrixXd::Constant(uc.rows(), uc.cols(), 1) - cp_dum;
-
/* Forces on panels */
- std::cout << "cp" << std::endl;
+ /*std::cout << "cp" << std::endl;
std::cout << cp << std::endl;
std::cout << "pan->s" << std::endl;
std::cout << pan->s << std::endl;
std::cout << "pan->nx" << std::endl;
std::cout << pan->nx << std::endl;
std::cout << "pan->Surface" << std::endl;
- std::cout << pan->Surface << std::endl;
+ std::cout << pan->Surface << std::endl; */
Eigen::MatrixXd fx = -cp.array()*pan->s.array()*pan->nx.array()/pan->Surface;
- Eigen::MatrixXd fy = -cp.array()*pan->s.array()*pan->nx.array()/pan->Surface;
- Eigen::MatrixXd fz = -cp.array()*pan->s.array()*pan->nx.array()/pan->Surface;
+ Eigen::MatrixXd fy = -cp.array()*pan->s.array()*pan->ny.array()/pan->Surface;
+ Eigen::MatrixXd fz = -cp.array()*pan->s.array()*pan->nz.array()/pan->Surface;
- std::cout << fx << std::endl;
+ /* std::cout << fx << std::endl;
std::cout << fy << std::endl;
- std::cout << fz << std::endl;
+ std::cout << fz << std::endl; */
double Fx = fx.sum();
//double Fy = fy.sum();
@@ -495,7 +510,18 @@ void SDSolver::Run()
CL(it) = Fz-Fx*fs_alpha;
CD(it) = Fz*fs_alpha+Fx;
}
- std::cout << " CL = " << CL(it) << std::endl;
+
+ double error = (CL(it)-CL(it-1))/CL(it-1);
+
+ std::cout << std::fixed << std::setprecision(5);
+ std::cout << std::setw(8) << it
+ << std::setw(12) << CL(it)
+ << std::setw(12) << CD(it)
+ << std::setw(15) << std::log10(error) << "\n";
+
+ if (error < tol){
+ break;
+ }
it++;
}
diff --git a/sdpm/src/pSDSolver.h b/sdpm/src/pSDSolver.h
index 9180461..8c6cef8 100644
--- a/sdpm/src/pSDSolver.h
+++ b/sdpm/src/pSDSolver.h
@@ -21,12 +21,13 @@ public:
Panels *pan;
double fs_alpha, fs_Mach, fs_airspeed;
unsigned long maxIt;
+double tol;
Eigen::VectorXd CL;
Eigen::VectorXd CD;
Eigen::VectorXd CY;
Eigen::MatrixXd cp;
-SDSolver(Panels *_pan, double _alpha, double _airspeed, unsigned long _maxIt);
+SDSolver(Panels *_pan, double _alpha, double _airspeed, double tol=1e-3, unsigned long _maxIt=100);
void Run();
diff --git a/sdpm/tests/panels3d.py b/sdpm/tests/panels3d.py
index 04bb3b1..b7c5056 100644
--- a/sdpm/tests/panels3d.py
+++ b/sdpm/tests/panels3d.py
@@ -8,19 +8,15 @@ def main():
nPanels = 2
airspeed = 50
- maxIt = 10
wing = sdpm.Panels(nPanels)
wing.CreateWing()
- solver = sdpm.SDSolver(wing, alpha, airspeed, maxIt)
+ solver = sdpm.SDSolver(wing, alpha, airspeed)
print('Computation started (pyCheck)')
solver.Run()
- plt.plot(solver.CL)
- plt.show()
-
print('Computation ended (pyCheck)')
if __name__ == "__main__":
diff --git a/vii/models/n0012_3_visc.geo b/vii/models/n0012_3_visc.geo
new file mode 100644
index 0000000..b449ddc
--- /dev/null
+++ b/vii/models/n0012_3_visc.geo
@@ -0,0 +1,397 @@
+/* Rectangular NACA0012 wing */
+// Initially generated by unsRgridWingGen.m
+// For gmsh 4, use line 334 instead of line 335 (Bezier <- BSpline)
+
+// Parameters
+// domain and mesh
+DefineConstant[ spn = { 2.0, Name "Wing span" } ];
+DefineConstant[ lgt = { 3.0, Name "Channel length" } ];
+DefineConstant[ wdt = { 3.0, Name "Channel width" } ];
+DefineConstant[ hgt = { 3.0, Name "Channel height" } ];
+DefineConstant[ msLe = { 0.05, Name "Leading edge mesh size" } ];
+DefineConstant[ msTe = { 0.05, Name "Trailing edge mesh size" } ];
+DefineConstant[ msF = { 1.0, Name "Farfield mesh size" } ];
+
+//// GEOMETRY
+
+
+/// Points
+
+// Airfoil 1: naca0012, 101 points
+Point(1) = {1.000000,0.000000,0.000000,msTe};
+Point(2) = {0.999013,0.000000,0.000143};
+Point(3) = {0.996057,0.000000,0.000572};
+Point(4) = {0.991144,0.000000,0.001280};
+Point(5) = {0.984292,0.000000,0.002260};
+Point(6) = {0.975528,0.000000,0.003501};
+Point(7) = {0.964888,0.000000,0.004990};
+Point(8) = {0.952414,0.000000,0.006710};
+Point(9) = {0.938153,0.000000,0.008643};
+Point(10) = {0.922164,0.000000,0.010770};
+Point(11) = {0.904508,0.000000,0.013071,msTe};
+Point(12) = {0.885257,0.000000,0.015523};
+Point(13) = {0.864484,0.000000,0.018106};
+Point(14) = {0.842274,0.000000,0.020795};
+Point(15) = {0.818712,0.000000,0.023569};
+Point(16) = {0.793893,0.000000,0.026405};
+Point(17) = {0.767913,0.000000,0.029279};
+Point(18) = {0.740877,0.000000,0.032168};
+Point(19) = {0.712890,0.000000,0.035048};
+Point(20) = {0.684062,0.000000,0.037896};
+Point(21) = {0.654508,0.000000,0.040686};
+Point(22) = {0.624345,0.000000,0.043394};
+Point(23) = {0.593691,0.000000,0.045992};
+Point(24) = {0.562667,0.000000,0.048455};
+Point(25) = {0.531395,0.000000,0.050754};
+Point(26) = {0.500000,0.000000,0.052862};
+Point(27) = {0.468605,0.000000,0.054749};
+Point(28) = {0.437333,0.000000,0.056390};
+Point(29) = {0.406309,0.000000,0.057755};
+Point(30) = {0.375655,0.000000,0.058819};
+Point(31) = {0.345492,0.000000,0.059557};
+Point(32) = {0.315938,0.000000,0.059947};
+Point(33) = {0.287110,0.000000,0.059971,msTe};
+Point(34) = {0.259123,0.000000,0.059614};
+Point(35) = {0.232087,0.000000,0.058863};
+Point(36) = {0.206107,0.000000,0.057712};
+Point(37) = {0.181288,0.000000,0.056159};
+Point(38) = {0.157726,0.000000,0.054206};
+Point(39) = {0.135516,0.000000,0.051862};
+Point(40) = {0.114743,0.000000,0.049138};
+Point(41) = {0.095492,0.000000,0.046049};
+Point(42) = {0.077836,0.000000,0.042615};
+Point(43) = {0.061847,0.000000,0.038859};
+Point(44) = {0.047586,0.000000,0.034803};
+Point(45) = {0.035112,0.000000,0.030473};
+Point(46) = {0.024472,0.000000,0.025893};
+Point(47) = {0.015708,0.000000,0.021088};
+Point(48) = {0.008856,0.000000,0.016078};
+Point(49) = {0.003943,0.000000,0.010884};
+Point(50) = {0.000987,0.000000,0.005521};
+Point(51) = {0.000000,0.000000,0.000000,msLe};
+Point(52) = {0.000987,0.000000,-0.005521};
+Point(53) = {0.003943,0.000000,-0.010884};
+Point(54) = {0.008856,0.000000,-0.016078};
+Point(55) = {0.015708,0.000000,-0.021088};
+Point(56) = {0.024472,0.000000,-0.025893};
+Point(57) = {0.035112,0.000000,-0.030473};
+Point(58) = {0.047586,0.000000,-0.034803};
+Point(59) = {0.061847,0.000000,-0.038859};
+Point(60) = {0.077836,0.000000,-0.042615};
+Point(61) = {0.095492,0.000000,-0.046049};
+Point(62) = {0.114743,0.000000,-0.049138};
+Point(63) = {0.135516,0.000000,-0.051862};
+Point(64) = {0.157726,0.000000,-0.054206};
+Point(65) = {0.181288,0.000000,-0.056159};
+Point(66) = {0.206107,0.000000,-0.057712};
+Point(67) = {0.232087,0.000000,-0.058863};
+Point(68) = {0.259123,0.000000,-0.059614};
+Point(69) = {0.287110,0.000000,-0.059971,msTe};
+Point(70) = {0.315938,0.000000,-0.059947};
+Point(71) = {0.345492,0.000000,-0.059557};
+Point(72) = {0.375655,0.000000,-0.058819};
+Point(73) = {0.406309,0.000000,-0.057755};
+Point(74) = {0.437333,0.000000,-0.056390};
+Point(75) = {0.468605,0.000000,-0.054749};
+Point(76) = {0.500000,0.000000,-0.052862};
+Point(77) = {0.531395,0.000000,-0.050754};
+Point(78) = {0.562667,0.000000,-0.048455};
+Point(79) = {0.593691,0.000000,-0.045992};
+Point(80) = {0.624345,0.000000,-0.043394};
+Point(81) = {0.654508,0.000000,-0.040686};
+Point(82) = {0.684062,0.000000,-0.037896};
+Point(83) = {0.712890,0.000000,-0.035048};
+Point(84) = {0.740877,0.000000,-0.032168};
+Point(85) = {0.767913,0.000000,-0.029279};
+Point(86) = {0.793893,0.000000,-0.026405};
+Point(87) = {0.818712,0.000000,-0.023569};
+Point(88) = {0.842274,0.000000,-0.020795};
+Point(89) = {0.864484,0.000000,-0.018106};
+Point(90) = {0.885257,0.000000,-0.015523};
+Point(91) = {0.904508,0.000000,-0.013071,msTe};
+Point(92) = {0.922164,0.000000,-0.010770};
+Point(93) = {0.938153,0.000000,-0.008643};
+Point(94) = {0.952414,0.000000,-0.006710};
+Point(95) = {0.964888,0.000000,-0.004990};
+Point(96) = {0.975528,0.000000,-0.003501};
+Point(97) = {0.984292,0.000000,-0.002260};
+Point(98) = {0.991144,0.000000,-0.001280};
+Point(99) = {0.996057,0.000000,-0.000572};
+Point(100) = {0.999013,0.000000,-0.000143,msTe};
+
+// Airfoil 2: naca0012, 101 points
+Point(102) = {1.000000,spn,0.000000,msTe};
+Point(103) = {0.999013,spn,0.000143};
+Point(104) = {0.996057,spn,0.000572};
+Point(105) = {0.991144,spn,0.001280};
+Point(106) = {0.984292,spn,0.002260};
+Point(107) = {0.975528,spn,0.003501};
+Point(108) = {0.964888,spn,0.004990};
+Point(109) = {0.952414,spn,0.006710};
+Point(110) = {0.938153,spn,0.008643};
+Point(111) = {0.922164,spn,0.010770};
+Point(112) = {0.904508,spn,0.013071,msTe};
+Point(113) = {0.885257,spn,0.015523};
+Point(114) = {0.864484,spn,0.018106};
+Point(115) = {0.842274,spn,0.020795};
+Point(116) = {0.818712,spn,0.023569};
+Point(117) = {0.793893,spn,0.026405};
+Point(118) = {0.767913,spn,0.029279};
+Point(119) = {0.740877,spn,0.032168};
+Point(120) = {0.712890,spn,0.035048};
+Point(121) = {0.684062,spn,0.037896};
+Point(122) = {0.654508,spn,0.040686};
+Point(123) = {0.624345,spn,0.043394};
+Point(124) = {0.593691,spn,0.045992};
+Point(125) = {0.562667,spn,0.048455};
+Point(126) = {0.531395,spn,0.050754};
+Point(127) = {0.500000,spn,0.052862};
+Point(128) = {0.468605,spn,0.054749};
+Point(129) = {0.437333,spn,0.056390};
+Point(130) = {0.406309,spn,0.057755};
+Point(131) = {0.375655,spn,0.058819};
+Point(132) = {0.345492,spn,0.059557};
+Point(133) = {0.315938,spn,0.059947};
+Point(134) = {0.287110,spn,0.059971,msTe};
+Point(135) = {0.259123,spn,0.059614};
+Point(136) = {0.232087,spn,0.058863};
+Point(137) = {0.206107,spn,0.057712};
+Point(138) = {0.181288,spn,0.056159};
+Point(139) = {0.157726,spn,0.054206};
+Point(140) = {0.135516,spn,0.051862};
+Point(141) = {0.114743,spn,0.049138};
+Point(142) = {0.095492,spn,0.046049};
+Point(143) = {0.077836,spn,0.042615};
+Point(144) = {0.061847,spn,0.038859};
+Point(145) = {0.047586,spn,0.034803};
+Point(146) = {0.035112,spn,0.030473};
+Point(147) = {0.024472,spn,0.025893};
+Point(148) = {0.015708,spn,0.021088};
+Point(149) = {0.008856,spn,0.016078};
+Point(150) = {0.003943,spn,0.010884};
+Point(151) = {0.000987,spn,0.005521};
+Point(152) = {0.000000,spn,0.000000,msLe};
+Point(153) = {0.000987,spn,-0.005521};
+Point(154) = {0.003943,spn,-0.010884};
+Point(155) = {0.008856,spn,-0.016078};
+Point(156) = {0.015708,spn,-0.021088};
+Point(157) = {0.024472,spn,-0.025893};
+Point(158) = {0.035112,spn,-0.030473};
+Point(159) = {0.047586,spn,-0.034803};
+Point(160) = {0.061847,spn,-0.038859};
+Point(161) = {0.077836,spn,-0.042615};
+Point(162) = {0.095492,spn,-0.046049};
+Point(163) = {0.114743,spn,-0.049138};
+Point(164) = {0.135516,spn,-0.051862};
+Point(165) = {0.157726,spn,-0.054206};
+Point(166) = {0.181288,spn,-0.056159};
+Point(167) = {0.206107,spn,-0.057712};
+Point(168) = {0.232087,spn,-0.058863};
+Point(169) = {0.259123,spn,-0.059614};
+Point(170) = {0.287110,spn,-0.059971,msTe};
+Point(171) = {0.315938,spn,-0.059947};
+Point(172) = {0.345492,spn,-0.059557};
+Point(173) = {0.375655,spn,-0.058819};
+Point(174) = {0.406309,spn,-0.057755};
+Point(175) = {0.437333,spn,-0.056390};
+Point(176) = {0.468605,spn,-0.054749};
+Point(177) = {0.500000,spn,-0.052862};
+Point(178) = {0.531395,spn,-0.050754};
+Point(179) = {0.562667,spn,-0.048455};
+Point(180) = {0.593691,spn,-0.045992};
+Point(181) = {0.624345,spn,-0.043394};
+Point(182) = {0.654508,spn,-0.040686};
+Point(183) = {0.684062,spn,-0.037896};
+Point(184) = {0.712890,spn,-0.035048};
+Point(185) = {0.740877,spn,-0.032168};
+Point(186) = {0.767913,spn,-0.029279};
+Point(187) = {0.793893,spn,-0.026405};
+Point(188) = {0.818712,spn,-0.023569};
+Point(189) = {0.842274,spn,-0.020795};
+Point(190) = {0.864484,spn,-0.018106};
+Point(191) = {0.885257,spn,-0.015523};
+Point(192) = {0.904508,spn,-0.013071,msTe};
+Point(193) = {0.922164,spn,-0.010770};
+Point(194) = {0.938153,spn,-0.008643};
+Point(195) = {0.952414,spn,-0.006710};
+Point(196) = {0.964888,spn,-0.004990};
+Point(197) = {0.975528,spn,-0.003501};
+Point(198) = {0.984292,spn,-0.002260};
+Point(199) = {0.991144,spn,-0.001280};
+Point(200) = {0.996057,spn,-0.000572};
+Point(201) = {0.999013,spn,-0.000143,msTe};
+
+// Tip:
+Point(5101) = {0.999013,spn,0.000000};
+Point(5102) = {0.996057,spn+0.000125,0.000000};
+Point(5103) = {0.991144,spn+0.000331,0.000000};
+Point(5104) = {0.984292,spn+0.000620,0.000000};
+Point(5105) = {0.975528,spn+0.000989,0.000000};
+Point(5106) = {0.964888,spn+0.001437,0.000000};
+Point(5107) = {0.952414,spn+0.001963,0.000000};
+Point(5108) = {0.938153,spn+0.002563,0.000000};
+Point(5109) = {0.922164,spn+0.003237,0.000000};
+Point(5110) = {0.904508,spn+0.003981,0.000000,msTe};
+Point(5111) = {0.885257,spn+0.004791,0.000000};
+Point(5112) = {0.864484,spn+0.005666,0.000000};
+Point(5113) = {0.842274,spn+0.006602,0.000000};
+Point(5114) = {0.818712,spn+0.007594,0.000000};
+Point(5115) = {0.793893,spn+0.008640,0.000000};
+Point(5116) = {0.767913,spn+0.009734,0.000000};
+Point(5117) = {0.740877,spn+0.010873,0.000000};
+Point(5118) = {0.712890,spn+0.012052,0.000000};
+Point(5119) = {0.684062,spn+0.013266,0.000000};
+Point(5120) = {0.654508,spn+0.014511,0.000000};
+Point(5121) = {0.624345,spn+0.015781,0.000000};
+Point(5122) = {0.593691,spn+0.017072,0.000000};
+Point(5123) = {0.562667,spn+0.018379,0.000000};
+Point(5124) = {0.531395,spn+0.019696,0.000000};
+Point(5125) = {0.500000,spn+0.021019,0.000000};
+Point(5126) = {0.468605,spn+0.022341,0.000000};
+Point(5127) = {0.437333,spn+0.023658,0.000000};
+Point(5128) = {0.406309,spn+0.024965,0.000000};
+Point(5129) = {0.375655,spn+0.026256,0.000000};
+Point(5130) = {0.345492,spn+0.027526,0.000000};
+Point(5131) = {0.315938,spn+0.028771,0.000000};
+Point(5132) = {0.287110,spn+0.029985,0.000000,msTe};
+Point(5133) = {0.019492,spn+0.041801,0.000000};
+
+// Dummy tip center:
+Point(5349) = {0.904508,spn,0.000000};
+Point(5350) = {0.287110,spn,0.000000};
+
+// Box
+Point(10001) = {1+lgt, 0, 0, msF};
+
+// Wake
+Point(10021) = {1+lgt, spn, 0, msF};
+
+/// Lines
+
+// Airfoil 1:
+Spline(1) = {1,2,3,4,5,6,7,8,9,10,11};
+Spline(2) = {11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33};
+Spline(3) = {33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51};
+Spline(4) = {51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69};
+Spline(5) = {69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91};
+Spline(6) = {91,92,93,94,95,96,97,98,99,100,1};
+
+// Airfoil 2:
+Spline(7) = {102,103,104,105,106,107,108,109,110,111,112};
+Spline(8) = {112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134};
+Spline(9) = {134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152};
+Spline(10) = {152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170};
+Spline(11) = {170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192};
+Spline(12) = {192,193,194,195,196,197,198,199,200,201,102};
+
+// Airfoil 1 to airfoil 2:
+Line(61) = {1,102};
+Line(62) = {11,112};
+Line(63) = {33,134};
+Line(64) = {51,152};
+Line(65) = {69,170};
+Line(66) = {91,192};
+
+// Tip:
+Spline(121) = {102,5101,5102,5103,5104,5105,5106,5107,5108,5109,5110};
+Spline(122) = {5110,5111,5112,5113,5114,5115,5116,5117,5118,5119,5120,5121,5122,5123,5124,5125,5126,5127,5128,5129,5130,5131,5132};
+If(GMSH_MAJOR_VERSION >= 4)
+ Bezier(123) = {5132,5133,152};
+Else
+ BSpline(123) = {5132,5133,152};
+EndIf
+Ellipse(124) = {112,5349,112,5110};
+Ellipse(125) = {134,5350,134,5132};
+Ellipse(126) = {170,5350,170,5132};
+Ellipse(127) = {192,5349,192,5110};
+
+// Box
+Line(221) = {1, 10001};
+Line(222) = {102, 10021};
+
+// Wake
+Line(241) = {10001, 10021};
+
+/// Line loops & Surfaces
+
+// Wing 1:
+Line Loop(1) = {1,62,-7,-61};
+Line Loop(2) = {2,63,-8,-62};
+Line Loop(3) = {3,64,-9,-63};
+Line Loop(4) = {4,65,-10,-64};
+Line Loop(5) = {5,66,-11,-65};
+Line Loop(6) = {6,61,-12,-66};
+Surface(1) = {-1};
+Surface(2) = {-2};
+Surface(3) = {-3};
+Surface(4) = {-4};
+Surface(5) = {-5};
+Surface(6) = {-6};
+
+// Wingtip:
+Line Loop(11) = {7,124,-121};
+Line Loop(12) = {8,125,-122,-124};
+Line Loop(13) = {9,-123,-125};
+Line Loop(14) = {10,126,123};
+Line Loop(15) = {11,127,122,-126};
+Line Loop(16) = {12,121,-127};
+Surface(11) = {-11};
+Surface(12) = {-12};
+Surface(13) = {-13};
+Surface(14) = {-14};
+Surface(15) = {-15};
+Surface(16) = {-16};
+
+// Wake
+Line Loop(51) = {221, 241, -222, -61};
+Surface(51) = {51};
+
+//// MESHING ALGORITHM
+//+
+nElmSpan = 50;
+Transfinite Curve {65, 63, 62, 66, 61, 64, 241} = nElmSpan Using Progression 1;
+nElmChrodwise_wake = 20;
+//+
+Transfinite Curve {222, 221} = nElmChrodwise_wake Using Progression 1.01;
+Transfinite Surface {51};
+
+nChordwise_midWing = 40;
+//+
+Transfinite Curve {8, 11, 5, 2, 122} = nChordwise_midWing Using Progression 1.01;
+//+
+nChordwise_Le = 40;
+Transfinite Curve {9, 9, 10, 3, 4} = nChordwise_Le Using Progression 0.95;
+
+//+
+nChordwise_Te = 20;
+Transfinite Curve {7, 121, 12, 1, 6} = nChordwise_Te Using Progression 1.01;
+
+Transfinite Surface {3};
+Transfinite Surface {4};
+Transfinite Surface {2};
+Transfinite Surface {5};
+Transfinite Surface {1};
+Transfinite Surface {6};
+
+Recombine Surface {3};
+Recombine Surface {4};
+Recombine Surface {2};
+Recombine Surface {5};
+Recombine Surface {1};
+Recombine Surface {6};
+Recombine Surface {51};
+
+
+//// PHYSICAL GROUPS
+
+// Trailing edge and wake tip
+Physical Line("wakeTip") = {222};
+Physical Line("te") = {61};
+Physical Line("teTip") = {61, 222};
+
+// Wing:
+Physical Surface("wing") = {1,2,3,11,12,13};
+Physical Surface("wing_") = {4,5,6,14,15,16};
+
+// Wake:
+Physical Surface("wake") = {51};
diff --git a/vii/pyVII/vInterpolator.py b/vii/pyVII/vInterpolator.py
index 66dee72..573cc10 100644
--- a/vii/pyVII/vInterpolator.py
+++ b/vii/pyVII/vInterpolator.py
@@ -1,367 +1,686 @@
-import numpy as np
-from matplotlib import pyplot as plt
-from fwk.coloring import ccolors
-
-import dart.pyVII.vUtils as blxU
-import sys
-#import dart.pyVII.vStructures as vStruct
-
-# 3D spline with Bézier curves
-from scipy.interpolate import bisplrep
-from scipy.interpolate import bisplev
-from scipy.interpolate import splrep
-from scipy.interpolate import splev
-
-from scipy.interpolate import interp1d
-
-from scipy.interpolate import RBFInterpolator
-
-import fwk
-
-
-class Interpolator:
- def __init__(self, _dartAPI, _blw, _cfg):
-
- self.iSolver = _dartAPI
- self.blw = _blw
-
- # Create data structures for interpolation
- self.tms = fwk.Timers()
- self.tms['GetWing'].start()
- self.__GetWing(_blw, _cfg)
- self.tms['GetWing'].stop()
-
- self.cfg = _cfg
- print(self.tms)
-
- self.stgPt=np.zeros(self.cfg['nSections'], dtype=int)
-
- self.xxPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(self.cfg['nSections'])]
- #self.uePrev = [[np.zeros(len(self.viscStruct[iSec][iReg][:,0])) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
- self.DeltaStarPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(self.cfg['nSections'])]
-
-
- def GetInvBC(self, vSolver, couplIter):
- self.tms['InvToVisc'].start()
-
- # Extract parameters from the inviscid solver
- # Particular care must be taken with space directions; In the inviscid, the airfoil is in the x,z plane.
- # While setting up the viscous solver, y and z must be inverted.
-
- # Wing
-
- U = [np.zeros((len(self.invStruct[i][:,0]), 3)) for i in range(len(self.invStruct))]
- M = [np.zeros(len(self.invStruct[i][:,0])) for i in range(len(self.invStruct))]
- Rho = [np.zeros(len(self.invStruct[i][:,0])) for i in range(len(self.invStruct))]
-
- # Extract variables.
-
- for iReg in range(len(self.invStruct)):
- for iNode, row in enumerate(self.invStruct[iReg][:,3]):
- row = int(row)
- U[iReg][iNode,:] = [self.iSolver.U[row][0], self.iSolver.U[row][1], self.iSolver.U[row][2]]
- M[iReg][iNode] = self.iSolver.M[row]
- Rho[iReg][iNode] = self.iSolver.rho[row]
-
- if self.cfg['nDim'] == 3:
- U[iReg][:,[1,2]] = U[iReg][:,[2,1]]
-
- # Find stagnation point
- nD = self.cfg['nDim']
- for iSec in range(len(self.viscStruct_tr)):
- Ux = [self.__RbfToSections(self.invStruct_tr[0][:,:nD], U[0][:,0], self.viscStruct_tr[iSec][0][:,:nD]), self.__RbfToSections(self.invStruct_tr[1][:, [0,2] if nD==3 else 0], U[1][:,0], self.viscStruct[iSec][1][:,[0,2] if nD==3 else 0])]
- Uy = [self.__RbfToSections(self.invStruct_tr[0][:,:nD], U[0][:,1], self.viscStruct_tr[iSec][0][:,:nD]), self.__RbfToSections(self.invStruct_tr[1][:, [0,2] if nD==3 else 0], U[1][:,1], self.viscStruct[iSec][1][:,[0,2] if nD==3 else 0])]
- Uz = [self.__RbfToSections(self.invStruct_tr[0][:,:nD], U[0][:,2], self.viscStruct_tr[iSec][0][:,:nD]), self.__RbfToSections(self.invStruct_tr[1][:, [0,2] if nD==3 else 0], U[1][:,2], self.viscStruct[iSec][1][:,[0,2] if nD==3 else 0])]
- Me = [self.__RbfToSections(self.invStruct_tr[0][:,:nD], M[0], self.viscStruct_tr[iSec][0][:,:nD]), self.__RbfToSections(self.invStruct_tr[1][:, [0,2] if nD==3 else 0], M[1], self.viscStruct[iSec][1][:,[0,2] if nD==3 else 0])]
- Rhoe = [self.__RbfToSections(self.invStruct_tr[0][:,:nD], Rho[0], self.viscStruct_tr[iSec][0][:,:nD]), self.__RbfToSections(self.invStruct_tr[1][:, [0,2] if nD==3 else 0], Rho[1], self.viscStruct[iSec][1][:,[0,2] if nD==3 else 0])]
-
- if couplIter == 0:
- self.stgPt[iSec] = np.argmin(Ux[0]**2+Uy[0]**2)
-
- xUp, xLw = self.__Remesh(self.viscStruct[iSec][0][:,0], self.stgPt[iSec])
- yUp, yLw = self.__Remesh(self.viscStruct[iSec][0][:,1], self.stgPt[iSec])
- zUp, zLw = self.__Remesh(self.viscStruct[iSec][0][:,2], self.stgPt[iSec])
-
- vSolver.SetGlobMesh(iSec, 0, xUp, yUp, zUp)
- vSolver.SetGlobMesh(iSec, 1, xLw, yLw, zLw)
- vSolver.SetGlobMesh(iSec, 2, self.viscStruct[iSec][1][:,0],
- self.viscStruct[iSec][1][:,1],
- self.viscStruct[iSec][1][:,2])
-
- self.DeltaStarPrev[iSec][0] = np.zeros(len(xUp))
- self.DeltaStarPrev[iSec][1] = np.zeros(len(xLw))
- self.DeltaStarPrev[iSec][2] = np.zeros(len(self.viscStruct[iSec][1][:,0]))
-
- self.xxPrev[iSec][0] = np.zeros(len(xUp))
- self.xxPrev[iSec][1] = np.zeros(len(xLw))
- self.xxPrev[iSec][2] = np.zeros(len(self.viscStruct[iSec][1][:,0]))
-
- UxUp, UxLw = self.__Remesh(Ux[0], self.stgPt[iSec])
- UyUp, UyLw = self.__Remesh(Uy[0], self.stgPt[iSec])
- UzUp, UzLw = self.__Remesh(Uz[0], self.stgPt[iSec])
-
- plt.plot(xUp, UxUp)
- plt.plot(xLw, UxLw)
- plt.show()
-
- vSolver.SetVelocities(iSec, 0, UxUp, UyUp, UzUp)
- vSolver.SetVelocities(iSec, 1, UxLw, UyLw, UzLw)
- vSolver.SetVelocities(iSec, 2, Ux[1], Uy[1], Uz[1])
-
- MeUp, MeLw = self.__Remesh(Me[0], self.stgPt[iSec])
-
- vSolver.SetMe(iSec, 0, MeUp)
- vSolver.SetMe(iSec, 1, MeLw)
- vSolver.SetMe(iSec, 2, Me[1])
-
- RhoUp, RhoLw = self.__Remesh(Rhoe[0], self.stgPt[iSec])
-
- vSolver.SetRhoe(iSec, 0, RhoUp)
- vSolver.SetRhoe(iSec, 1, RhoLw)
- vSolver.SetRhoe(iSec, 2, Rhoe[1])
-
- for iReg in range(3):
- vSolver.SetDeltaStarExt(iSec, iReg, self.DeltaStarPrev[iSec][iReg])
- vSolver.SetxxExt(iSec, iReg, self.xxPrev[iSec][iReg])
-
-
-
- def SetBlowingVel(self, vSolver, couplIter):
-
- for iSec in range(len(self.viscStruct)):
- for iReg in range(3):
- self.DeltaStarPrev[iSec][iReg] = vSolver.ExtractDeltaStar(iSec, iReg)
- self.xxPrev[iSec][iReg] = vSolver.Extractxx(iSec, iReg)
-
- x = [np.zeros((0,3)) for _ in range(2)]
- blw = [np.zeros(0) for _ in range(2)]
-
- for iSec in range(len(self.viscStruct)):
- for iReg in range(2):
- if iReg == 0:
- x[iReg] = np.vstack((x[iReg], self.viscCg_tr[iSec][iReg][:,:3]))
- blwUp = vSolver.ExtractBlowingVel(iSec, 0)
- blwLw = vSolver.ExtractBlowingVel(iSec, 1)
- blw[iReg] = np.concatenate((blw[iReg], blwUp[::-1] , blwLw))
- elif iReg == 1:
- x[iReg] = np.vstack((x[iReg], self.viscCg_tr[iSec][iReg][:,:3]))
- blw[iReg] = np.concatenate((blw[iReg], vSolver.ExtractBlowingVel(iSec, 2)))
-
- self.tms['Interpolation'].start()
- nD = self.cfg['nDim']
- mappedBlw = [self.__RbfToSections(x[0], blw[0], self.invCg_tr[0][:,:3]), self.__RbfToSections(x[1][:, [0,2] if nD==3 else 0], blw[1], self.invCg_tr[1][:, [0,2] if nD==3 else 0])]
-
- self.tms['Interpolation'].stop()
- for iElm in range(len(self.invCg_tr[0][:,0])):
- self.blw[0].setU(iElm, mappedBlw[0][iElm])
- for iElm in range(len(self.invCg_tr[1][:,0])):
- self.blw[1].setU(iElm, mappedBlw[1][iElm])
-
- def RunSolver(self):
- self.iSolver.run()
-
- def ResetSolver(self):
- self.iSolver.reset()
-
- def GetCp(self,bnd):
- import dart.utils as floU
- Cp = floU.extract(bnd, self.iSolver.Cp)
- return Cp
-
- def GetAlpha(self):
- return self.iSolver.pbl.alpha
-
- def GetMinf(self):
- return self.iSolver.pbl.M_inf
-
- def GetCl(self):
- return self.iSolver.Cl
-
- def GetCm(self):
- return self.iSolver.Cm
-
- def __Remesh(self, vec, stgPt):
- xUp = vec[:stgPt+1]
- xUp = xUp[::-1]
-
- xLw = vec[stgPt:]
- return xUp, xLw
-
- def __RbfToSections(self, x, var, xs):
- if np.all(var == var[0]):
- varOut = RBFInterpolator(x, var, neighbors=1, kernel='linear', smoothing=0.0, degree=0)(xs)
- else:
- varOut = RBFInterpolator(x, var, neighbors=self.cfg['neighbors'], kernel=self.cfg['rbftype'], smoothing=self.cfg['smoothing'], degree=self.cfg['degree'])(xs)
- return varOut
-
-
- def __GetWing(self, blw, config):
-
- viscStruct = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(2)] for _ in range(config['nSections'])]
- viscStruct_tr = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(2)] for _ in range(config['nSections'])]
-
-
- invStruct = [np.zeros((len(blw[i].nodes), 4)) for i in range(len(blw))]
- invStruct_tr = [np.zeros((len(blw[i].nodes), 4)) for i in range(len(blw))]
-
- # Wing
- upperNodes = np.zeros((0,4))
- lowerNodes = np.zeros((0,4))
-
- self.tms['ElmLoop'].start()
- # Compute reference normal
- nVec = [np.zeros(0) for _ in range(blw[1].tag.elems[0].nodes.size())]
- for iNode, n in enumerate(blw[1].tag.elems[0].nodes):
- nVec[iNode] = np.array([n.pos[0], n.pos[1], n.pos[2]])
- normalRef = np.cross(nVec[1] - nVec[0], nVec[2] - nVec[0])
- normalRef/=np.linalg.norm(normalRef)
-
- for e in blw[0].tag.elems:
- nVec = [np.zeros(0) for _ in range(e.nodes.size())]
- for iNode, n in enumerate(e.nodes):
- nVec[iNode] = np.array([n.pos[0], n.pos[1], n.pos[2]])
- normal = np.cross(nVec[1] - nVec[0], nVec[2] - nVec[0])
- normal/=np.linalg.norm(normalRef)
- if np.dot(normal, normalRef) >= 0:
- for n in e.nodes:
- if n.row not in upperNodes[:,3] and n.row not in lowerNodes[:,3]:
- upperNodes = np.vstack((upperNodes, [n.pos[0], n.pos[1], n.pos[2], n.row]))
- if np.dot(normal, normalRef) < 0:
- for n in e.nodes:
- if n.row not in lowerNodes[:,3] and n.row not in upperNodes[:,3]:
- lowerNodes = np.vstack((lowerNodes, [n.pos[0], n.pos[1], n.pos[2], n.row]))
-
- self.tms['ElmLoop'].stop()
-
- dummyLower = lowerNodes.copy()
- dummyLower[:,0] *= -1
-
- invStruct_tr[0] = np.vstack((upperNodes, dummyLower))
- invStruct[0] = np.vstack((upperNodes, lowerNodes))
-
- # Remove duplicates
- #invStruct_tr[0] = np.unique(invStruct_tr[0], axis=0)
-
- # Wake part
- for iNode, n in enumerate(blw[1].nodes):
- invStruct_tr[1][iNode, :] = [n.pos[0], n.pos[1], n.pos[2], n.row]
- invStruct[1] = invStruct_tr[1].copy()
-
- # Create spanwise distribution
- if config['nDim'] == 2:
- spanDistri = np.zeros(1)
- spanDistri[0] = invStruct[0][0,2] # Span direction is the z direction
-
- wingSpan = 0.
- fInterp = interp1d(invStruct_tr[0][:,0], invStruct_tr[0][:,1])
-
- elif config['nDim'] == 3:
- invStruct[0][:, [1,2]] = invStruct[0][:, [2, 1]]
- invStruct[1][:, [1,2]] = invStruct[1][:, [2, 1]]
- invStruct_tr[0][:, [1,2]] = invStruct_tr[0][:, [2, 1]]
- invStruct_tr[1][:, [1,2]] = invStruct_tr[1][:, [2, 1]]
- spanDistri = np.linspace(min(invStruct[0][:,2]), max(invStruct[0][:,2]), config['nSections'])
- spanDistri[0] +=0.1
- spanDistri[-1] -=0.1
-
- # Spline wing
- # Point in the wing tip region are not accounted for during Bezier curve spline
-
- m = len(invStruct_tr[0][:,0])
- _s = (m-np.sqrt(2*m) + m+np.sqrt(2*m))/2
- tck_wing = bisplrep(invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],0], invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],2], invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],1], kx=5, ky=5, s=1e-6, quiet=1)
-
- # Create sections
- self.tms['iSecLoop'].start()
-
- for iSec, z in enumerate(spanDistri):
- portion = invStruct[0][abs(invStruct[0][:,2] - z) <= 0.01 * config['span'],:]
-
- chordLength = max(portion[:, 0]) - min(portion[:, 0])
- xDistri = min(portion[:,0]) + 1/2*(1 + np.cos(np.linspace(0, np.pi, config['nPoints'][0])))*chordLength
-
- viscStruct[iSec][0][:,0] = np.concatenate((xDistri[:-1], xDistri[::-1]), axis=None)
- viscStruct[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
-
- viscStruct_tr[iSec][0][:,0] = np.concatenate((xDistri[:-1], -xDistri[::-1]), axis=None).copy()
- viscStruct_tr[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
-
- # Interpolation in the transformed space.
- if config['nDim'] == 3:
-
- for iPoint in range(len(viscStruct_tr[iSec][0][:,0])):
- viscStruct_tr[iSec][0][iPoint, 1] = bisplev(viscStruct_tr[iSec][0][iPoint,0], viscStruct_tr[iSec][0][iPoint,2], tck_wing)
-
-
- elif config['nDim'] == 2:
- viscStruct_tr[iSec][0][:,1] = fInterp(viscStruct_tr[iSec][0][:,0])
-
- #viscStruct_tr[iSec][0][0,1] = 0.
- #viscStruct_tr[iSec][0][-1,1] = 0.
- viscStruct[iSec][0][:,1] = viscStruct_tr[iSec][0][:,1].copy()
-
- # Wake
- wakeLength = np.max(invStruct[1][:,0]) - np.max(viscStruct[iSec][0][:,0])
- viscStruct_tr[iSec][1][:,0] = np.max(viscStruct[iSec][0][:,0]) + (np.max(viscStruct[iSec][0][:,0]) + np.cos(np.linspace(np.pi, np.pi/2, config['nPoints'][1]))) * wakeLength
- viscStruct_tr[iSec][1][:,2] = z*np.ones(len(viscStruct[iSec][1][:,0]))
- viscStruct_tr[iSec][1][:,1] = interp1d(invStruct_tr[1][:,0], invStruct_tr[1][:,1])(viscStruct_tr[iSec][1][:,0])
-
- viscStruct[iSec][1] = viscStruct_tr[iSec][1].copy()
-
- # Cg positions
- self.tms['iSecLoop'].stop()
-
- """fig = plt.figure()
- ax = fig.add_subplot(projection='3d')
- ax.scatter(invStruct_tr[0][:,0], invStruct_tr[0][:,2], invStruct_tr[0][:,1], s=0.3)
- ax.scatter(invStruct_tr[1][:,0], invStruct_tr[1][:,2], invStruct_tr[1][:,1], s=0.3)
- for sec in viscStruct_tr:
- ax.scatter(sec[0][:,0], sec[0][:,2], sec[0][:,1], color = 'red')
- ax.scatter(sec[1][:,0], sec[1][:,2], sec[1][:,1], color = 'red')
- ax.set_zlim([-0.5, 0.5])
- ax.set_xlabel('x')
- ax.set_ylabel('y')
- ax.set_zlabel('z')
- plt.show()"""
-
- viscCg_tr = [[np.zeros((len(viscStruct[iSec][iReg][:,0]) - 1, 3)) for iReg in range(2)] for iSec in range(len(viscStruct))]
- invCg_tr = [np.zeros((len(blw[iReg].tag.elems), 3)) for iReg in range(2)]
- self.tms['CgLoop'].start()
-
- # Viscous
- for iSec in range(len(viscStruct)):
- for iReg in range(len(viscStruct[iSec])):
- for iElm in range(len(viscStruct[iSec][iReg][:,0]) - 1):
- for iDir in range(3):
- viscCg_tr[iSec][iReg][iElm, iDir] = viscStruct_tr[iSec][iReg][iElm, iDir] + viscStruct_tr[iSec][iReg][iElm + 1, iDir]
- viscCg_tr[iSec][iReg] /= 2
-
- # Inviscid
- for iReg in range(len(blw)):
- for iElm, e in enumerate(blw[iReg].tag.elems):
- for n in e.nodes:
- if n.row in upperNodes[:,3] or iReg == 1:
- invCg_tr[iReg][iElm, 0] +=n.pos[0]
- elif n.row in lowerNodes[:,3]:
- invCg_tr[iReg][iElm, 0] += -n.pos[0]
-
- invCg_tr[iReg][iElm, 1] +=n.pos[1]
- invCg_tr[iReg][iElm, 2] +=n.pos[2]
- invCg_tr[iReg][iElm, :] /= e.nodes.size()
-
- if config['nDim'] == 3:
- invCg_tr[iReg][:,[1,2]] = invCg_tr[iReg][:,[2,1]]
- self.tms['CgLoop'].stop()
-
- self.invStruct = invStruct
- self.invStruct_tr = invStruct_tr
- self.viscStruct = viscStruct
- self.viscStruct_tr = viscStruct_tr
-
- self.invCg_tr = invCg_tr
- self.viscCg_tr = viscCg_tr
- np.set_printoptions(threshold=sys.maxsize)
- for sec in self.viscStruct:
- print('next')
- sec[0][:,[1,2]] = sec[0][:,[2,1]]
- print(sec[0])
+import numpy as np
+from matplotlib import pyplot as plt
+from fwk.coloring import ccolors
+
+# 3D spline with Bézier curves
+from scipy.interpolate import bisplrep
+from scipy.interpolate import bisplev
+from scipy.interpolate import splrep
+from scipy.interpolate import splev
+
+from scipy.interpolate import interp1d
+
+from scipy.interpolate import RBFInterpolator
+
+import fwk
+
+
+class Interpolator:
+ def __init__(self, _dartAPI, _blw, _imsh, _vmsh, _cfg):
+
+ self.cfg = _cfg
+
+ self.tms = fwk.Timers()
+
+ self.tms['Mesh sort'].start()
+ self.idata, self.icg = self.GetMesh(_imsh)
+ self.vdata, self.vcg = self.GetMesh(_vmsh)
+ self.sortedRows, self.sortedElems = self.SortMesh(self.vdata, self.vcg)
+
+ self.tms['Mesh sort'].stop()
+ print(self.tms)
+
+ fig = plt.figure()
+ ax = fig.add_subplot(projection='3d')
+ ax.scatter(self.vdata[0][:,0], self.vdata[0][:,1], self.vdata[0][:,2], s=0.3)
+ ax.scatter(self.vcg[0][:,0], self.vcg[0][:,1], self.vcg[0][:,2], s=0.3)
+ ax.scatter(self.vdata[1][:,0], self.vdata[1][:,1], self.vdata[1][:,2], s=0.3)
+ ax.scatter(self.vcg[1][:,0], self.vcg[1][:,1], self.vcg[1][:,2], s=0.3)
+ ax.set_xlabel('x')
+ ax.set_ylabel('y')
+ ax.set_zlabel('z')
+ #plt.gca().invert_yaxis()
+ plt.axis('off')
+ plt.show()
+
+ fig = plt.figure()
+ ax = fig.add_subplot(projection='3d')
+ ax.scatter(self.idata[0][:,0], self.idata[0][:,1], self.idata[0][:,2], s=0.3)
+ ax.scatter(self.icg[0][:,0], self.icg[0][:,1], self.icg[0][:,2], s=0.3)
+ ax.scatter(self.idata[1][:,0], self.idata[1][:,1], self.idata[1][:,2], s=0.3)
+ ax.scatter(self.icg[1][:,0], self.icg[1][:,1], self.icg[1][:,2], s=0.3)
+ ax.set_xlabel('x')
+ ax.set_ylabel('y')
+ ax.set_zlabel('z')
+ #plt.gca().invert_yaxis()
+ plt.axis('off')
+ plt.show()
+
+ self.iSolver = _dartAPI
+ self.blw = _blw
+
+ # Create data structures for interpolation
+ self.tms = fwk.Timers()
+ self.tms['GetWing'].start()
+ """if _cfg['nDim'] == 2:
+ self.__GetAirfoil(_blw, _cfg)
+ elif _cfg['nDim'] == 3:
+ self.__GetWing(_blw, _cfg)"""
+ self.tms['GetWing'].stop()
+
+ # self.cfg = _cfg
+ print(self.tms)
+
+ self.stgPt=np.zeros(self.cfg['nSections'], dtype=int)
+
+ self.xxPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(self.cfg['nSections'])]
+ #self.uePrev = [[np.zeros(len(self.viscStruct[iSec][iReg][:,0])) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
+ self.DeltaStarPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(self.cfg['nSections'])]
+
+
+ def GetInvBC(self, vSolver, couplIter):
+ self.tms['InvToVisc'].start()
+
+ # Extract parameters from the inviscid solver
+ # Particular care must be taken with space directions; In the inviscid, the airfoil is in the x,z plane.
+ # While setting up the viscous solver, y and z must be inverted.
+
+ # Wing
+
+ U = [np.zeros((len(self.idata[i][:,0]), 3)) for i in range(len(self.idata))]
+ M = [np.zeros(len(self.idata[i][:,0])) for i in range(len(self.idata))]
+ Rho = [np.zeros(len(self.idata[i][:,0])) for i in range(len(self.idata))]
+
+ # Extract variables.
+
+ for iReg in range(len(self.idata)):
+ for iNode, row in enumerate(self.idata[iReg][:,3]):
+ row = int(row)
+ U[iReg][iNode,:] = [self.iSolver.U[row][0], self.iSolver.U[row][1], self.iSolver.U[row][2]]
+ M[iReg][iNode] = self.iSolver.M[row]
+ Rho[iReg][iNode] = self.iSolver.rho[row]
+
+ """if self.cfg['nDim'] == 3:
+ U[iReg][:,[1,2]] = U[iReg][:,[2,1]]"""
+
+ # Find stagnation point
+ nD = self.cfg['nDim']
+ Ux = [self.__RbfToSections(self.idata[0][:,:nD], U[0][:,0], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, [0,2] if nD==3 else 0], U[1][:,0], self.vdata[1][:,[0,2] if nD==3 else 0])]
+ Uy = [self.__RbfToSections(self.idata[0][:,:nD], U[0][:,1], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, [0,2] if nD==3 else 0], U[1][:,1], self.vdata[1][:,[0,2] if nD==3 else 0])]
+ Uz = [self.__RbfToSections(self.idata[0][:,:nD], U[0][:,2], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, [0,2] if nD==3 else 0], U[1][:,2], self.vdata[1][:,[0,2] if nD==3 else 0])]
+ Me = [self.__RbfToSections(self.idata[0][:,:nD], M[0], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, [0,2] if nD==3 else 0], M[1], self.vdata[1][:,[0,2] if nD==3 else 0])]
+ Rhoe = [self.__RbfToSections(self.idata[0][:,:nD], Rho[0], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, [0,2] if nD==3 else 0], Rho[1], self.vdata[1][:,[0,2] if nD==3 else 0])]
+
+ print('yeeeeah')
+ quit()
+ if couplIter == 0:
+ self.stgPt[iSec] = np.argmin(Ux[0]**2+Uy[0]**2)
+
+ xUp, xLw = self.__Remesh(self.viscStruct[iSec][0][:,0], self.stgPt[iSec])
+ yUp, yLw = self.__Remesh(self.viscStruct[iSec][0][:,1], self.stgPt[iSec])
+ zUp, zLw = self.__Remesh(self.viscStruct[iSec][0][:,2], self.stgPt[iSec])
+
+ vSolver.SetGlobMesh(iSec, 0, xUp, yUp, zUp)
+ vSolver.SetGlobMesh(iSec, 1, xLw, yLw, zLw)
+ vSolver.SetGlobMesh(iSec, 2, self.viscStruct[iSec][1][:,0],
+ self.viscStruct[iSec][1][:,1],
+ self.viscStruct[iSec][1][:,2])
+
+ self.DeltaStarPrev[iSec][0] = np.zeros(len(xUp))
+ self.DeltaStarPrev[iSec][1] = np.zeros(len(xLw))
+ self.DeltaStarPrev[iSec][2] = np.zeros(len(self.viscStruct[iSec][1][:,0]))
+
+ self.xxPrev[iSec][0] = np.zeros(len(xUp))
+ self.xxPrev[iSec][1] = np.zeros(len(xLw))
+ self.xxPrev[iSec][2] = np.zeros(len(self.viscStruct[iSec][1][:,0]))
+
+ UxUp, UxLw = self.__Remesh(Ux[0], self.stgPt[iSec])
+ UyUp, UyLw = self.__Remesh(Uy[0], self.stgPt[iSec])
+ UzUp, UzLw = self.__Remesh(Uz[0], self.stgPt[iSec])
+
+ plt.plot(xUp, UxUp)
+ plt.plot(xLw, UxLw)
+ plt.xlabel('$x/c$')
+ plt.ylabel('$U_x$')
+ plt.show()
+
+ vSolver.SetVelocities(iSec, 0, UxUp, UyUp, UzUp)
+ vSolver.SetVelocities(iSec, 1, UxLw, UyLw, UzLw)
+ vSolver.SetVelocities(iSec, 2, Ux[1], Uy[1], Uz[1])
+
+ MeUp, MeLw = self.__Remesh(Me[0], self.stgPt[iSec])
+
+ vSolver.SetMe(iSec, 0, MeUp)
+ vSolver.SetMe(iSec, 1, MeLw)
+ vSolver.SetMe(iSec, 2, Me[1])
+
+ RhoUp, RhoLw = self.__Remesh(Rhoe[0], self.stgPt[iSec])
+
+ vSolver.SetRhoe(iSec, 0, RhoUp)
+ vSolver.SetRhoe(iSec, 1, RhoLw)
+ vSolver.SetRhoe(iSec, 2, Rhoe[1])
+
+ for iReg in range(3):
+ vSolver.SetDeltaStarExt(iSec, iReg, self.DeltaStarPrev[iSec][iReg])
+ vSolver.SetxxExt(iSec, iReg, self.xxPrev[iSec][iReg])
+
+
+
+ def SetBlowingVel(self, vSolver, couplIter):
+
+ for iSec in range(len(self.viscStruct)):
+ for iReg in range(3):
+ self.DeltaStarPrev[iSec][iReg] = vSolver.ExtractDeltaStar(iSec, iReg)
+ self.xxPrev[iSec][iReg] = vSolver.Extractxx(iSec, iReg)
+
+ x = [np.zeros((0,3)) for _ in range(2)]
+ blw = [np.zeros(0) for _ in range(2)]
+
+ for iSec in range(len(self.viscStruct)):
+ for iReg in range(2):
+ if iReg == 0:
+ x[iReg] = np.vstack((x[iReg], self.viscCg_tr[iSec][iReg][:,:3]))
+ blwUp = vSolver.ExtractBlowingVel(iSec, 0)
+ blwLw = vSolver.ExtractBlowingVel(iSec, 1)
+ blw[iReg] = np.concatenate((blw[iReg], blwUp[::-1] , blwLw))
+ elif iReg == 1:
+ x[iReg] = np.vstack((x[iReg], self.viscCg_tr[iSec][iReg][:,:3]))
+ blw[iReg] = np.concatenate((blw[iReg], vSolver.ExtractBlowingVel(iSec, 2)))
+
+ fig = plt.figure()
+ ax = fig.add_subplot(projection='3d')
+ ax.scatter(x[0][:,0], x[0][:,2], blw[0], color = 'red')
+ """for sec in viscStruct_tr:
+ ax.scatter(sec[0][:,0], sec[0][:,2], sec[0][:,1], color = 'red')
+ #ax.scatter(sec[1][:,0], sec[1][:,2], sec[1][:,1], color = 'red')"""
+ #ax.set_zlim([-0.5, 0.5])
+ ax.set_xlabel('x')
+ ax.set_ylabel('y')
+ ax.set_zlabel('Blowing velocity')
+ #plt.axis('off')
+ plt.show()
+ self.tms['Interpolation'].start()
+ nD = self.cfg['nDim']
+ #wignggg = self.__RbfToSections(x[0], blw[0], self.invCg_tr[0][:,:3])
+ wakke = self.__RbfToSections(x[1][:, [0,2] if nD==3 else 0], blw[1], self.invCg_tr[1][:, [0,2] if nD==3 else 0])
+ mappedBlw = [self.__RbfToSections(x[0], blw[0], self.invCg_tr[0][:,:3]), self.__RbfToSections(x[1][:, [0,2] if nD==3 else 0], blw[1], self.invCg_tr[1][:, [0,2] if nD==3 else 0])]
+
+ self.tms['Interpolation'].stop()
+ for iElm in range(len(self.invCg_tr[0][:,0])):
+ self.blw[0].setU(iElm, mappedBlw[0][iElm])
+ for iElm in range(len(self.invCg_tr[1][:,0])):
+ self.blw[1].setU(iElm, mappedBlw[1][iElm])
+
+ def RunSolver(self):
+ self.iSolver.run()
+
+ def ResetSolver(self):
+ self.iSolver.reset()
+
+ def GetCp(self,bnd):
+ import dart.utils as floU
+ Cp = floU.extract(bnd, self.iSolver.Cp)
+ return Cp
+
+ def GetAlpha(self):
+ return self.iSolver.pbl.alpha
+
+ def GetMinf(self):
+ return self.iSolver.pbl.M_inf
+
+ def GetCl(self):
+ return self.iSolver.Cl
+
+ def GetCm(self):
+ return self.iSolver.Cm
+
+ def __Remesh(self, vec, stgPt):
+ xUp = vec[:stgPt+1]
+ xUp = xUp[::-1]
+
+ xLw = vec[stgPt:]
+ return xUp, xLw
+
+ def __RbfToSections(self, x, var, xs):
+ if np.all(var == var[0]):
+ varOut = RBFInterpolator(x, var, neighbors=1, kernel='linear', smoothing=0.0, degree=0)(xs)
+ else:
+ varOut = RBFInterpolator(x, var, neighbors=self.cfg['neighbors'], kernel=self.cfg['rbftype'], smoothing=self.cfg['smoothing'], degree=self.cfg['degree'])(xs)
+ return varOut
+
+
+ def __GetWing(self, blw, config):
+
+ viscStruct = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(2)] for _ in range(config['nSections'])]
+ viscStruct_tr = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(2)] for _ in range(config['nSections'])]
+
+
+ invStruct = [np.zeros((len(blw[i].nodes), 4)) for i in range(len(blw))]
+ invStruct_tr = [np.zeros((len(blw[i].nodes), 4)) for i in range(len(blw))]
+
+ # Wing
+ upperNodes = np.zeros((0,4))
+ lowerNodes = np.zeros((0,4))
+
+ self.tms['ElmLoop'].start()
+ # Compute reference normal
+ nVec = [np.zeros(0) for _ in range(blw[1].tag.elems[0].nodes.size())]
+ for iNode, n in enumerate(blw[1].tag.elems[0].nodes):
+ nVec[iNode] = np.array([n.pos[0], n.pos[1], n.pos[2]])
+ normalRef = np.cross(nVec[1] - nVec[0], nVec[2] - nVec[0])
+ normalRef/=np.linalg.norm(normalRef)
+
+ for e in blw[0].tag.elems:
+ nVec = [np.zeros(0) for _ in range(e.nodes.size())]
+ for iNode, n in enumerate(e.nodes):
+ nVec[iNode] = np.array([n.pos[0], n.pos[1], n.pos[2]])
+ normal = np.cross(nVec[1] - nVec[0], nVec[2] - nVec[0])
+ normal/=np.linalg.norm(normalRef)
+ if np.dot(normal, normalRef) >= 0:
+ for n in e.nodes:
+ if n.row not in upperNodes[:,3] and n.row not in lowerNodes[:,3]:
+ upperNodes = np.vstack((upperNodes, [n.pos[0], n.pos[1], n.pos[2], n.row]))
+ if np.dot(normal, normalRef) < 0:
+ for n in e.nodes:
+ if n.row not in lowerNodes[:,3] and n.row not in upperNodes[:,3]:
+ lowerNodes = np.vstack((lowerNodes, [n.pos[0], n.pos[1], n.pos[2], n.row]))
+
+ self.tms['ElmLoop'].stop()
+
+ dummyLower = lowerNodes.copy()
+ dummyLower[:,0] *= -1
+
+ invStruct_tr[0] = np.vstack((upperNodes, dummyLower))
+ invStruct[0] = np.vstack((upperNodes, lowerNodes))
+
+ # Remove duplicates
+ #invStruct_tr[0] = np.unique(invStruct_tr[0], axis=0)
+
+ # Wake part
+ for iNode, n in enumerate(blw[1].nodes):
+ invStruct_tr[1][iNode, :] = [n.pos[0], n.pos[1], n.pos[2], n.row]
+ invStruct[1] = invStruct_tr[1].copy()
+
+ # Create spanwise distribution
+ if config['nDim'] == 2:
+ spanDistri = np.zeros(1)
+ spanDistri[0] = invStruct[0][0,2] # Span direction is the z direction
+
+ wingSpan = 0.
+ fInterp = interp1d(invStruct_tr[0][:,0], invStruct_tr[0][:,1])
+
+ elif config['nDim'] == 3:
+ invStruct[0][:, [1,2]] = invStruct[0][:, [2, 1]]
+ invStruct[1][:, [1,2]] = invStruct[1][:, [2, 1]]
+ invStruct_tr[0][:, [1,2]] = invStruct_tr[0][:, [2, 1]]
+ invStruct_tr[1][:, [1,2]] = invStruct_tr[1][:, [2, 1]]
+ spanDistri = np.linspace(min(invStruct[0][:,2]), max(invStruct[0][:,2]), config['nSections'])
+ spanDistri[0] +=0.1
+ spanDistri[-1] -=0.1
+
+ # Spline wing
+ # Point in the wing tip region are not accounted for during Bezier curve spline
+
+ m = len(invStruct_tr[0][:,0])
+ _s = (m-np.sqrt(2*m) + m+np.sqrt(2*m))/2
+ tck_wing = bisplrep(invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],0], invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],2], invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],1], kx=5, ky=5, s=1e-6, quiet=1)
+
+ # Create sections
+ self.tms['iSecLoop'].start()
+
+ for iSec, z in enumerate(spanDistri):
+ portion = invStruct[0][abs(invStruct[0][:,2] - z) <= 0.01 * config['span'],:]
+
+ chordLength = max(portion[:, 0]) - min(portion[:, 0])
+ xDistri = min(portion[:,0]) + 1/2*(1 + np.cos(np.linspace(0, np.pi, config['nPoints'][0])))*chordLength
+
+ viscStruct[iSec][0][:,0] = np.concatenate((xDistri[:-1], xDistri[::-1]), axis=None)
+ viscStruct[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
+
+ viscStruct_tr[iSec][0][:,0] = np.concatenate((xDistri[:-1], -xDistri[::-1]), axis=None).copy()
+ viscStruct_tr[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
+
+ # Interpolation in the transformed space.
+ if config['nDim'] == 3:
+
+ for iPoint in range(len(viscStruct_tr[iSec][0][:,0])):
+ viscStruct_tr[iSec][0][iPoint, 1] = bisplev(viscStruct_tr[iSec][0][iPoint,0], viscStruct_tr[iSec][0][iPoint,2], tck_wing)
+
+
+ elif config['nDim'] == 2:
+ viscStruct_tr[iSec][0][:,1] = fInterp(viscStruct_tr[iSec][0][:,0])
+
+ #viscStruct_tr[iSec][0][0,1] = 0.
+ #viscStruct_tr[iSec][0][-1,1] = 0.
+ viscStruct[iSec][0][:,1] = viscStruct_tr[iSec][0][:,1].copy()
+
+ # Wake
+ wakeLength = np.max(invStruct[1][:,0]) - np.max(viscStruct[iSec][0][:,0])
+ viscStruct_tr[iSec][1][:,0] = np.max(viscStruct[iSec][0][:,0]) + (np.max(viscStruct[iSec][0][:,0]) + np.cos(np.linspace(np.pi, np.pi/2, config['nPoints'][1]))) * wakeLength
+ viscStruct_tr[iSec][1][:,2] = z*np.ones(len(viscStruct[iSec][1][:,0]))
+ viscStruct_tr[iSec][1][:,1] = interp1d(invStruct_tr[1][:,0], invStruct_tr[1][:,1])(viscStruct_tr[iSec][1][:,0])
+
+ viscStruct[iSec][1] = viscStruct_tr[iSec][1].copy()
+
+ # Cg positions
+ self.tms['iSecLoop'].stop()
+
+ viscCg_tr = [[np.zeros((len(viscStruct[iSec][iReg][:,0]) - 1, 3)) for iReg in range(2)] for iSec in range(len(viscStruct))]
+ invCg_tr = [np.zeros((len(blw[iReg].tag.elems), 3)) for iReg in range(2)]
+ self.tms['CgLoop'].start()
+
+ # Viscous
+ for iSec in range(len(viscStruct)):
+ for iReg in range(len(viscStruct[iSec])):
+ for iElm in range(len(viscStruct[iSec][iReg][:,0]) - 1):
+ for iDir in range(3):
+ viscCg_tr[iSec][iReg][iElm, iDir] = viscStruct_tr[iSec][iReg][iElm, iDir] + viscStruct_tr[iSec][iReg][iElm + 1, iDir]
+ viscCg_tr[iSec][iReg] /= 2
+
+ # Inviscid
+ for iReg in range(len(blw)):
+ for iElm, e in enumerate(blw[iReg].tag.elems):
+ for n in e.nodes:
+ if n.row in upperNodes[:,3] or iReg == 1:
+ invCg_tr[iReg][iElm, 0] +=n.pos[0]
+ elif n.row in lowerNodes[:,3]:
+ invCg_tr[iReg][iElm, 0] += -n.pos[0]
+
+ invCg_tr[iReg][iElm, 1] +=n.pos[1]
+ invCg_tr[iReg][iElm, 2] +=n.pos[2]
+ invCg_tr[iReg][iElm, :] /= e.nodes.size()
+
+ if config['nDim'] == 3:
+ invCg_tr[iReg][:,[1,2]] = invCg_tr[iReg][:,[2,1]]
+ self.tms['CgLoop'].stop()
+
+ self.invStruct = invStruct
+ self.invStruct_tr = invStruct_tr
+ self.viscStruct = viscStruct
+ self.viscStruct_tr = viscStruct_tr
+
+ self.invCg_tr = invCg_tr
+ self.viscCg_tr = viscCg_tr
+ """np.set_printoptions(threshold=sys.maxsize)
+ for sec in self.viscStruct:
+ print('next')
+ sec[0][:,[1,2]] = sec[0][:,[2,1]]
+ print(sec[0])"""
+
+ fig = plt.figure()
+ ax = fig.add_subplot(projection='3d')
+ ax.scatter(invStruct_tr[0][:,0], invStruct_tr[0][:,2], invStruct_tr[0][:,1], s=0.3)
+ ax.scatter(invStruct_tr[0][:,0], -invStruct_tr[0][:,2], invStruct_tr[0][:,1], s=0.3, color='grey')
+ #ax.scatter(invStruct[1][:,0], invStruct[1][:,2], invStruct[1][:,1], s=0.3)
+ for sec in viscStruct_tr:
+ ax.scatter(sec[0][:,0], sec[0][:,2], sec[0][:,1], color = 'red')
+ ax.scatter(sec[1][:,0], sec[1][:,2], sec[1][:,1], color = 'red')
+ ax.set_zlim([-0.5, 0.5])
+ ax.set_xlabel('x')
+ ax.set_ylabel('y')
+ ax.set_zlabel('z')
+ #plt.gca().invert_yaxis()
+ plt.axis('off')
+ plt.show()
+
+
+ def __GetAirfoil(self, blw, config):
+
+ if config['nSections'] != 1:
+ raise RuntimeError('Cannot create more than 1 section on 2D geometry. Specified:', config['nSections'])
+
+ viscStruct = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(2)] for _ in range(config['nSections'])]
+ viscStruct_tr = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(1)] for _ in range(config['nSections'])]
+
+ # Create node list
+ nodesCoord = np.zeros((0, 4))
+ for n in blw[0].nodes:
+ nodesCoord = np.vstack((nodesCoord, [n.pos[0], n.pos[1], n.pos[2], n.row]))
+
+ spanDistri = np.zeros(1)
+ spanDistri[0] = nodesCoord[0,2] # Span direction is the z direction
+
+ wingSpan = max(nodesCoord[:,2]) - min(nodesCoord[:,2])
+
+ #upperNodes = nodesCoord[nodesCoord[:,2]>0]
+
+ """nVec = [np.zeros(0) for _ in range(blw[1].tag.elems[0].nodes.size())]
+ nVec = blw[1].tag.elems[0]
+ normalRef = np.cross(nVec[1] - nVec[0], nVec[2] - nVec[0])
+ normalRef/=np.linalg.norm(normalRef)
+
+ upperNodes = np.zeros((0,4))
+ lowerNodes = np.zeros((0,4))
+ for e in blw[0].tag.elems:
+ # Compute normal
+ cVec = [e.nodes[1].pos[0] - e.nodes[0].pos[0], -(e.nodes[1].pos[1] - e.nodes[0].pos[1]), 0]
+ print(cVec)
+ if np.dot(cVec, normalRef) >=0:
+ for n in e.nodes:
+ if n.row not in upperNodes[:,3] and n.row not in lowerNodes[:,3]:
+ upperNodes = np.vstack((upperNodes, [n.pos[0], n.pos[1], n.pos[2], n.row]))
+
+ else:
+ for n in e.nodes:
+ if n.row not in upperNodes[:,3] and n.row not in lowerNodes[:,3]:
+ upperNodes = np.vstack((upperNodes, [n.pos[0], n.pos[1], n.pos[2], n.row]))
+
+ plt.plot(upperNodes[:,0], upperNodes[:,1])
+ plt.plot(lowerNodes[:,0], lowerNodes[:,1])
+ plt.show()"""
+ for iSec, z in enumerate(spanDistri):
+
+ portion = nodesCoord[abs(nodesCoord[:,2] - z) <= 0.01 * wingSpan,:]
+
+ le_coord = portion[portion[:,0] == np.min((portion[:,0])), :][0]
+ te_coord = portion[portion[:,0] == np.max((portion[:,0])), :][0]
+
+ # Find te elements
+ te_elems = []
+ for e in blw[0].tag.elems:
+ for n in e.nodes:
+ if n.pos[0] == te_coord[0]:
+ te_elems.append(e)
+ break
+
+ # Cg of te elements
+ cg_teElems = np.zeros((len(te_elems), 3))
+ for iElm, eTe in enumerate(te_elems):
+ for n in eTe.nodes:
+ cg_teElems[iElm, 0] += n.pos[0]
+ cg_teElems[iElm, 1] += n.pos[1]
+ cg_teElems[iElm, 2] += n.pos[2]
+
+ # Sort upper side in 0 and lower side in 1
+ if cg_teElems[0, 1] < cg_teElems[1, 1]:
+ save = te_elems[0].copy()
+ te_elems[0] = te_elems[1]
+ te_elems[1] = save
+
+ # Find te_nodes (sorted [upper, lower])
+ te_nodes = [None, None]
+ for i in range(len(te_elems)):
+ for n in te_elems[i].nodes:
+ if n.pos[0] == te_coord[0]:
+ te_nodes[i] = n
+ break
+
+ # Find le_node
+ for n in blw[0].nodes:
+ if n.pos[0] == le_coord[0]:
+ le_node = n
+ break
+
+
+ meany = (le_coord[1] + te_coord[1]) / 2
+
+ transformedCoord = np.zeros((0, 4))
+ nodesList = []
+
+ for iNode, n in enumerate(blw[0].nodes):
+ if n.pos[0] != le_coord[0] and n.pos[0] != te_coord[0]:
+ if n.pos[1] > meany:
+ transformedCoord = np.vstack((transformedCoord, [n.pos[0], n.pos[1], n.pos[2], n.row]))
+ nodesList.append(n)
+
+ elif n.pos[1] < meany:
+ transformedCoord = np.vstack((transformedCoord, [-n.pos[0], n.pos[1], n.pos[2], n.row]))
+ nodesList.append(n)
+
+ transformedCoord = np.vstack((transformedCoord, le_coord))
+ nodesList.append(le_node)
+ transformedCoord = np.vstack((transformedCoord, [te_nodes[0].pos[0], te_nodes[0].pos[1], te_nodes[0].pos[2], te_nodes[0].row]))
+ nodesList.append(te_nodes[0])
+ transformedCoord = np.vstack((transformedCoord, [-te_nodes[1].pos[0], te_nodes[1].pos[1], te_nodes[1].pos[2], te_nodes[1].row]))
+ nodesList.append(te_nodes[1])
+
+ wakeCoord = np.zeros((len(blw[1].nodes), 4))
+ for iNode, n in enumerate(blw[1].nodes):
+ wakeCoord[iNode, 0] = n.pos[0]
+ wakeCoord[iNode, 1] = n.pos[1]
+ wakeCoord[iNode, 2] = n.pos[2]
+ wakeCoord[iNode, 3] = n.row
+
+ chordLength = max(portion[:, 0]) - min(portion[:, 0])
+ xDistri = le_coord[0] + 1/2*(1 + np.cos(np.linspace(0, np.pi, config['nPoints'][0])))*chordLength
+
+ viscStruct[iSec][0][:,0] = np.concatenate((xDistri[:-1], xDistri[::-1]), axis=None)
+ viscStruct[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
+
+ viscStruct_tr[iSec][0][:,0] = np.concatenate((xDistri[:-1], -xDistri[::-1]), axis=None).copy()
+ viscStruct_tr[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
+
+ # Interpolation in the transformed space.
+ splineTr = transformedCoord.copy()
+ splineTr = splineTr[splineTr[:, 0].argsort()]
+ spl = splrep(splineTr[:,0], splineTr[:,1], k=5)
+ viscStruct_tr[iSec][0][:,1] = splev(viscStruct_tr[iSec][0][:,0], spl)
+ #viscStruct_tr[iSec][0][:,1] = interp1d(transformedCoord[:,0], transformedCoord[:,1])(viscStruct_tr[iSec][0][:,0])
+ viscStruct[iSec][0][:,1] = viscStruct_tr[iSec][0][:,1].copy()
+
+ # Wake
+ wakeLength = np.max(wakeCoord[:,0]) - np.max(viscStruct[iSec][0][:,0])
+ viscStruct[iSec][1][:,0] = np.max(viscStruct[iSec][0][:,0]) + (np.max(viscStruct[iSec][0][:,0]) + np.cos(np.linspace(np.pi, np.pi/2, config['nPoints'][1]))) * wakeLength
+ viscStruct[iSec][1][:,2] = z*np.ones(len(viscStruct[iSec][1][:,0]))
+ viscStruct[iSec][1][:,1] = interp1d(wakeCoord[:,0], wakeCoord[:,1])(viscStruct[iSec][1][:,0])
+
+ self.invStruct = [np.zeros((0, 4)) for _ in range(2)]
+ self.invStruct_tr = [np.zeros((0, 4)) for _ in range(2)]
+ self.invStruct[0] = abs(transformedCoord.copy())
+ self.invStruct_tr[0] = transformedCoord.copy()
+ self.invStruct[1] = abs(wakeCoord.copy())
+ self.invStruct_tr[1] = wakeCoord.copy()
+
+ print(np.shape(self.invStruct[0]))
+
+
+ self.viscStruct = viscStruct
+ self.viscStruct_tr = viscStruct_tr
+ #plt.plot(nodesCoord[:,0], nodesCoord[:,1], 'o', markersize=8, color='green')
+ plt.plot(viscStruct[0][0][:,0], viscStruct[0][0][:,1], 'x', color='blue')
+ plt.axis('off')
+ plt.show()
+ return nodesList, transformedCoord, wakeCoord, viscStruct, viscStruct_tr
+
+ def GetMesh(self, msh):
+
+ data = [np.zeros((0,4)) for _ in range(2)]
+ cg = [np.zeros((0,4)) for _ in range(2)]
+
+ alreadysaid = ['']
+ for e in msh.elems:
+ if e.ptag.name not in alreadysaid:
+ print(e.ptag.name)
+ alreadysaid.append(e.ptag.name)
+ print('')
+
+ for e in msh.elems:
+ # Compute cg position
+ cg_pt = np.zeros(3)
+ for n in e.nodes:
+ cg_pt += [n.pos[0], n.pos[1], n.pos[2]]
+ cg_pt/=e.nodes.size()
+ cg_pt = np.hstack((cg_pt, e.no))
+
+ if e.ptag.name == "wake":
+ for nw in e.nodes:
+ if nw.row not in data[1][:,3]:
+ data[1] = np.vstack((data[1], [nw.pos[0], nw.pos[1], nw.pos[2], nw.row]))
+ cg[1] = np.vstack((cg[1], cg_pt))
+ elif e.ptag.name == 'wing':
+ for nw in e.nodes:
+ if nw.row not in data[0][:,3]:
+ data[0] = np.vstack((data[0], [nw.pos[0], nw.pos[1], nw.pos[2], nw.row]))
+ cg[0] = np.vstack((cg[0], cg_pt))
+ elif e.ptag.name == 'wing_':
+ for nw in e.nodes:
+ if nw.row not in data[0][:,3]:
+ data[0] = np.vstack((data[0], [-nw.pos[0], nw.pos[1], nw.pos[2], nw.row]))
+ cg_pt[0] = -cg_pt[0]
+ cg[0] = np.vstack((cg[0], cg_pt))
+
+ return data, cg
+
+ def SortMesh(self, data, cg):
+
+ data = data.copy()
+
+ data[0] = data[0][data[0][:,1]<=1, :]
+
+ cg = data.copy()
+
+ cg[0] = cg[0][cg[0][:,1]<=1, :]
+
+ sections = data[0][:,1]
+ for i in range(len(sections)):
+ sections[i] = round(sections[i], 3)
+ sections = np.unique(sections)
+
+ nChord = [len(data[i][abs(data[i][:,1]-sections[0])<1e-3, 0]) for i in range(len(data))]
+
+ sortedRows = [np.zeros((nChord[i], 0)) for i in range(len(data))]
+
+ for iReg, reg in enumerate(data):
+
+ for iy, y in enumerate(sections):
+
+ # Find nearest value in the mesh
+
+ idx = (np.abs(reg[:,1]-y).argmin())
+ sections[iy] = reg[idx,1]
+
+ sort = reg[abs(reg[:,1]-sections[iy])<1e-3, :]
+
+ sortedRows[iReg] = np.column_stack((sortedRows[iReg], sort[:,3]))
+
+ sections_cg = cg[0][:,1]
+ for i in range(len(sections_cg)):
+ sections_cg[i] = round(sections_cg[i], 3)
+ sections_cg = np.unique(sections)
+
+ nChord_cg = [len(cg[i][abs(cg[i][:,1]-sections[0])<1e-3, 0]) for i in range(len(cg))]
+
+ sortedElems = [np.zeros((nChord_cg[i], 0)) for i in range(len(cg))]
+
+ for iReg, reg in enumerate(cg):
+
+ for iy, y in enumerate(sections_cg):
+
+ # Find nearest value in the mesh
+
+ idx = (np.abs(reg[:,1]-y).argmin())
+ sections_cg[iy] = reg[idx,1]
+
+ sort = reg[abs(reg[:,1]-sections_cg[iy])<1e-3, :]
+
+ sortedElems[iReg] = np.column_stack((sortedElems[iReg], sort[:,3]))
+
+ return sortedRows, sortedElems
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/vii/pyVII/vUtils.py b/vii/pyVII/vUtils.py
index 16b8868..0a49e49 100644
--- a/vii/pyVII/vUtils.py
+++ b/vii/pyVII/vUtils.py
@@ -23,6 +23,30 @@ def initBL(Re, Minf, CFL0, nSections, verbose=1):
return solver
+def mesh(file, pars):
+ """Initialize mesh and mesh writer
+
+ Parameters
+ ----------
+ dim : int (2 or 3)
+ number of dimensions
+ file : str
+ file contaning mesh (w/o extention)
+ pars : dict
+ parameters for mesh
+ bnds : array of str
+ names of crack (wake) boundaries physical tag
+ wk : str (optional)
+ name of crack (wake) physical tag
+ wktp : str (optional)
+ name of crack (wake) tip physical tag
+ """
+ import tbox.gmsh as tmsh
+ # load the mesh
+ msh = tmsh.MeshLoader(file,__file__).execute(**pars)
+ return msh
+
+
def GetSolution(iSec, vSolver):
if iSec<0:
raise RuntimeError("dart::viscU Invalid section number", iSec)
diff --git a/vii/src/wBLRegion.cpp b/vii/src/wBLRegion.cpp
index deb501f..9f18080 100644
--- a/vii/src/wBLRegion.cpp
+++ b/vii/src/wBLRegion.cpp
@@ -94,16 +94,17 @@ void BLRegion::SetStagBC(double Re)
DeltaStar[0] = U[0] * U[1];
- std::vector<double> ClosParam = closSolver->LaminarClosures(U[0], U[1], Ret[0], invBnd->GetMe(0), name);
-
- Hstar[0] = ClosParam[0];
- Hstar2[0] = ClosParam[1];
- Hk[0] = ClosParam[2];
- Cf[0] = ClosParam[3];
- Cd[0] = ClosParam[4];
- Delta[0] = ClosParam[5];
- Cteq[0] = ClosParam[6];
- Us[0] = ClosParam[7];
+ std::vector<double> lamParam(8, 0.);
+ closSolver->LaminarClosures(lamParam, U[0], U[1], Ret[0], invBnd->GetMe(0), name);
+
+ Hstar[0] = lamParam[0];
+ Hstar2[0] = lamParam[1];
+ Hk[0] = lamParam[2];
+ Cf[0] = lamParam[3];
+ Cd[0] = lamParam[4];
+ Delta[0] = lamParam[5];
+ Cteq[0] = lamParam[6];
+ Us[0] = lamParam[7];
} // End SetStagBC
@@ -133,16 +134,17 @@ void BLRegion::SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<d
/* Laminar closures. */
- std::vector<double> ClosParam = closSolver->LaminarClosures(U[0], U[1], Ret[0], invBnd->GetMe(0), name);
-
- Hstar[0] = ClosParam[0];
- Hstar2[0] = ClosParam[1];
- Hk[0] = ClosParam[2];
- Cf[0] = ClosParam[3];
- Cd[0] = ClosParam[4];
- Delta[0] = ClosParam[5];
- Cteq[0] = ClosParam[6];
- Us[0] = ClosParam[7];
+ std::vector<double> lamParam(8, 0.);
+ closSolver->LaminarClosures(lamParam, U[0], U[1], Ret[0], invBnd->GetMe(0), name);
+
+ Hstar[0] = lamParam[0];
+ Hstar2[0] = lamParam[1];
+ Hk[0] = lamParam[2];
+ Cf[0] = lamParam[3];
+ Cd[0] = lamParam[4];
+ Delta[0] = lamParam[5];
+ Cteq[0] = lamParam[6];
+ Us[0] = lamParam[7];
for (size_t k = 0; k < mesh->GetnMarkers(); ++k)
{
diff --git a/vii/src/wClosures.cpp b/vii/src/wClosures.cpp
index 49c1aa8..0c0434a 100644
--- a/vii/src/wClosures.cpp
+++ b/vii/src/wClosures.cpp
@@ -12,7 +12,7 @@ Closures::~Closures(){}
/**
* @brief Laminar closure relations
*/
-std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret, double Me, std::string name) const
+void Closures::LaminarClosures(std::vector<double> &closureVars, double theta, double H, double Ret, const double Me, const std::string name) const
{
H = std::max(H, 1.00005);
@@ -90,29 +90,15 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
double Us = 0.;
double Cteq = 0.;
- std::vector<double> ClosureVars(8, 0.);
-
- ClosureVars[0] = Hstar;
- ClosureVars[1] = Hstar2;
- ClosureVars[2] = Hk;
- ClosureVars[3] = Cf;
- ClosureVars[4] = Cd;
- ClosureVars[5] = delta;
- ClosureVars[6] = Cteq;
- ClosureVars[7] = Us;
-
- return ClosureVars;
+ closureVars = {Hstar, Hstar2, Hk, Cf, Cd, delta, Cteq, Us};
}
/**
* @brief Turbulent closure relations
*/
-std::vector<double> Closures::TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name) const
+void Closures::TurbulentClosures(std::vector<double> &closureVars, double theta, double H, double Ct, double Ret, const double Me, const std::string name) const
{
- /*if (std::isnan(Ret))
- {
- std::cout << "Ret is nan. H:" << H << std::endl;
- }*/
+
H = std::max(H, 1.00005);
Ct = std::min(Ct, 0.3);
//Ct = std::max(std::min(Ct, 0.30), 0.0000001);
@@ -154,6 +140,7 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
{
Hstar = (Hk-H0)*(Hk-H0)*(0.007*std::log(Ret)/((Hk-H0+4/std::log(Ret))*(Hk-H0+4/std::log(Ret))) + 0.015/Hk) + 1.5 + 4/Ret;
}
+
/* Whitfield's minor additional compressibility correction */
Hstar = (Hstar + 0.028*Me*Me)/(1+0.014*Me*Me);
@@ -223,20 +210,92 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
std::cout << "Negative sqrt encountered " << std::endl;
//throw;
}
- //double Cteq = std::sqrt(n*Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H));
+
+ /* Dissipation coefficient */
Cd = n*(Cdw+ Cdd + Cdl);
- // Ueq = 1.25/Hk*(Cf/2-((Hk-1)/(6.432*Hk))**2)
- std::vector<double> ClosureVars(8, 0.);
+ closureVars = {Hstar, Hstar2, Hk, Cf, Cd, delta, Cteq, Us};
+
+}
+
+void Closures::TurbulentClosures(double &closureVars, double theta, double H, double Ct, double Ret, const double Me, const std::string name) const
+{
- ClosureVars[0] = Hstar;
- ClosureVars[1] = Hstar2;
- ClosureVars[2] = Hk;
- ClosureVars[3] = Cf;
- ClosureVars[4] = Cd;
- ClosureVars[5] = delta;
- ClosureVars[6] = Cteq;
- ClosureVars[7] = Us;
+ H = std::max(H, 1.00005);
+ Ct = std::min(Ct, 0.3);
- return ClosureVars;
+ double Hk = (H - 0.29*Me*Me)/(1+0.113*Me*Me);
+ if (name == "wake")
+ {
+ Hk = std::max(Hk, 1.00005);
+ }
+ else
+ {
+ Hk = std::max(Hk, 1.05000);
+ }
+
+ double Hstar2 = ((0.064/(Hk-0.8))+0.251)*Me*Me;
+ double gamma = 1.4 - 1.;
+ double Fc = std::sqrt(1+0.5*gamma*Me*Me); /* .2 can be used as well replacing .5*gamma */
+
+ double H0 = 0.;
+ if (Ret <= 400)
+ {
+ H0 = 4.0;
+ }
+ else
+ {
+ H0 = 3 + (400/Ret);
+ }
+ if (Ret <= 200)
+ {
+ Ret = 200;
+ }
+
+ double Hstar = 0.;
+ if (Hk <= H0)
+ {
+ Hstar = ((0.5-4/Ret)*((H0-Hk)/(H0-1))*((H0-Hk)/(H0-1)))*(1.5/(Hk+0.5))+1.5+4/Ret;
+ }
+ else
+ {
+ Hstar = (Hk-H0)*(Hk-H0)*(0.007*std::log(Ret)/((Hk-H0+4/std::log(Ret))*(Hk-H0+4/std::log(Ret))) + 0.015/Hk) + 1.5 + 4/Ret;
+ }
+
+ /* Whitfield's minor additional compressibility correction */
+ Hstar = (Hstar + 0.028*Me*Me)/(1+0.014*Me*Me);
+
+ /* Equivalent normalized wall slip velocity */
+ double Us = (Hstar/2)*(1-4*(Hk-1)/(3*H));
+
+ double Ctcon = 0.5/(6.7*6.7*0.75);
+
+ double Hkc = 0.;
+ double Cteq = 0.;
+
+ if (name == "wake")
+ {
+ if (Us > 0.99995)
+ {
+ Us = 0.99995;
+ }
+ Hkc = Hk - 1;
+ Cteq = std::sqrt(4*Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H)); // Here it is Cteq^0.5
+ }
+ else
+ {
+ if (Us > 0.95)
+ {
+ Us = 0.98;
+ }
+ Hkc = std::max(Hk-1-18/Ret, 0.01);
+ Cteq = std::sqrt(Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H)); // Here it is Cteq^0.5
+ }
+ if (Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H)<0)
+ {
+ std::cout << "Negative sqrt encountered " << std::endl;
+ }
+
+ closureVars = Cteq;
+
}
\ No newline at end of file
diff --git a/vii/src/wClosures.h b/vii/src/wClosures.h
index 23f92d9..a786a79 100644
--- a/vii/src/wClosures.h
+++ b/vii/src/wClosures.h
@@ -13,8 +13,9 @@ class DART_API Closures
public:
Closures();
~Closures();
- std::vector<double> LaminarClosures(double theta, double H, double Ret, double Me, std::string name) const;
- std::vector<double> TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name) const;
+ void LaminarClosures(std::vector<double> &closureVars, double theta, double H, double Ret, const double Me, const std::string name) const;
+ void TurbulentClosures(std::vector<double> &closureVars, double theta, double H, double Ct, double Ret, double Me, const std::string name) const;
+ void TurbulentClosures(double &closureVars, double theta, double H, double Ct, double Ret, const double Me, const std::string name) const;
};
} // namespace vii
#endif //WClosures_H
\ No newline at end of file
diff --git a/vii/src/wTimeSolver.cpp b/vii/src/wTimeSolver.cpp
index be6e0e7..8c8ea39 100644
--- a/vii/src/wTimeSolver.cpp
+++ b/vii/src/wTimeSolver.cpp
@@ -194,7 +194,8 @@ void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double>
if (bl->Regime[iPoint] == 0)
{
- std::vector<double> lamParam = bl->closSolver->LaminarClosures(bl->U[iPoint * nVar + 0], bl->U[iPoint * nVar + 1], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
+ std::vector<double> lamParam(8, 0.);
+ bl->closSolver->LaminarClosures(lamParam, bl->U[iPoint * nVar + 0], bl->U[iPoint * nVar + 1], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
bl->Hstar[iPoint] = lamParam[0];
bl->Hstar2[iPoint] = lamParam[1];
bl->Hk[iPoint] = lamParam[2];
@@ -206,7 +207,8 @@ void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double>
} // End if
else if (bl->Regime[iPoint] == 1)
{
- std::vector<double> turbParam = bl->closSolver->TurbulentClosures(bl->U[iPoint * nVar + 0], bl->U[iPoint * nVar + 1], bl->U[iPoint * nVar + 4], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
+ std::vector<double> turbParam(8, 0.);
+ bl->closSolver->TurbulentClosures(turbParam, bl->U[iPoint * nVar + 0], bl->U[iPoint * nVar + 1], bl->U[iPoint * nVar + 4], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
bl->Hstar[iPoint] = turbParam[0];
bl->Hstar2[iPoint] = turbParam[1];
bl->Hk[iPoint] = turbParam[2];
diff --git a/vii/src/wViscFluxes.cpp b/vii/src/wViscFluxes.cpp
index 1b667ba..672aa99 100644
--- a/vii/src/wViscFluxes.cpp
+++ b/vii/src/wViscFluxes.cpp
@@ -83,7 +83,8 @@ VectorXd ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<double> u)
if (bl->Regime[iPoint] == 0)
{
/* Laminar closure relations */
- std::vector<double> lamParam = bl->closSolver->LaminarClosures(u[0], u[1], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
+ std::vector<double> lamParam(8, 0.);
+ bl->closSolver->LaminarClosures(lamParam, u[0], u[1], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
bl->Hstar[iPoint] = lamParam[0];
bl->Hstar2[iPoint] = lamParam[1];
bl->Hk[iPoint] = lamParam[2];
@@ -97,7 +98,8 @@ VectorXd ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<double> u)
else if (bl->Regime[iPoint] == 1)
{
/* Turbulent closure relations */
- std::vector<double> turbParam = bl->closSolver->TurbulentClosures(u[0], u[1], u[4], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
+ std::vector<double> turbParam(8, 0.);
+ bl->closSolver->TurbulentClosures(turbParam, u[0], u[1], u[4], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
bl->Hstar[iPoint] = turbParam[0];
bl->Hstar2[iPoint] = turbParam[1];
bl->Hk[iPoint] = turbParam[2];
diff --git a/vii/src/wViscSolver.cpp b/vii/src/wViscSolver.cpp
index c7dcb05..3e364fc 100644
--- a/vii/src/wViscSolver.cpp
+++ b/vii/src/wViscSolver.cpp
@@ -327,8 +327,9 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
double avLam = 1. - avTurb;
/* Impose boundary condition */
- std::vector<double> turbParam1 = bl->closSolver->TurbulentClosures(bl->U[(iPoint - 1) * nVar + 0], bl->U[(iPoint - 1) * nVar + 1], 0.03, bl->Ret[iPoint - 1], bl->invBnd->GetMe(iPoint - 1), bl->name);
- bl->Cteq[iPoint - 1] = turbParam1[6];
+ double Cteq_trans;
+ bl->closSolver->TurbulentClosures(Cteq_trans, bl->U[(iPoint - 1) * nVar + 0], bl->U[(iPoint - 1) * nVar + 1], 0.03, bl->Ret[iPoint - 1], bl->invBnd->GetMe(iPoint - 1), bl->name);
+ bl->Cteq[iPoint - 1] = Cteq_trans;
bl->U[(iPoint - 1) * nVar + 4] = 0.7 * bl->Cteq[iPoint - 1];
/* Avoid starting with ill conditioned IC for turbulent computation. (Regime was switched above) */
@@ -353,7 +354,8 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
bl->U[iPoint * nVar + 4] = avTurb * bl->U[(iPoint)*nVar + 4]; /* Ct */
/* Recompute closures. (The turbulent BC @ iPoint - 1 does not influence laminar closure relations). */
- std::vector<double> lamParam = bl->closSolver->LaminarClosures(bl->U[(iPoint - 1) * nVar + 0], bl->U[(iPoint - 1) * nVar + 1], bl->Ret[iPoint - 1], bl->invBnd->GetMe(iPoint - 1), bl->name);
+ std::vector<double> lamParam(8, 0.);
+ bl->closSolver->LaminarClosures(lamParam, bl->U[(iPoint - 1) * nVar + 0], bl->U[(iPoint - 1) * nVar + 1], bl->Ret[iPoint - 1], bl->invBnd->GetMe(iPoint - 1), bl->name);
bl->Hstar[iPoint - 1] = lamParam[0];
bl->Hstar2[iPoint - 1] = lamParam[1];
bl->Hk[iPoint - 1] = lamParam[2];
@@ -363,7 +365,8 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
bl->Cteq[iPoint - 1] = lamParam[6];
bl->Us[iPoint - 1] = lamParam[7];
- std::vector<double> turbParam = bl->closSolver->TurbulentClosures(bl->U[(iPoint)*nVar + 0], bl->U[(iPoint)*nVar + 1], bl->U[(iPoint)*nVar + 4], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
+ std::vector<double> turbParam(8, 0.);
+ bl->closSolver->TurbulentClosures(turbParam, bl->U[(iPoint)*nVar + 0], bl->U[(iPoint)*nVar + 1], bl->U[(iPoint)*nVar + 4], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
bl->Hstar[iPoint] = turbParam[0];
bl->Hstar2[iPoint] = turbParam[1];
bl->Hk[iPoint] = turbParam[2];
diff --git a/vii/tests/bli2.py b/vii/tests/bli2.py
index cf3b476..c46a10a 100644
--- a/vii/tests/bli2.py
+++ b/vii/tests/bli2.py
@@ -57,7 +57,7 @@ def main():
# define flow variables
Re = 1e7
- alpha = 2.*math.pi/180
+ alpha = 5.*math.pi/180
M_inf = 0.
CFL0 = 1
diff --git a/vii/tests/bli3.py b/vii/tests/bli3.py
index 718ba48..29c99f8 100644
--- a/vii/tests/bli3.py
+++ b/vii/tests/bli3.py
@@ -1,117 +1,120 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-## Compute lifting (linear or nonlinear) potential flow around a NACA 0012 rectangular wing
-# Adrien Crovato
-#
-# Test the nonlinear shock-capturing capability and the Kutta condition
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-import math
-import dart.utils as floU
-import dart.default as floD
-
-import dart.pyVII.vUtils as viscU
-import dart.pyVII.vCoupler2 as viscC
-import dart.pyVII.vInterpolator as vInterpol
-
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-def main():
- # timer
- tms = fwk.Timers()
- tms['total'].start()
-
- # define flow variables
- """Re = 11.72e6
- alpha = 3.06*math.pi/180
- M_inf = 0.8395"""
- Re = 1e7
- alpha = 0.*math.pi/180
- M_inf = 0.
- dim = 3
- CFL0 = 1
-
- nSections = 5
-
- # define dimension and mesh size
- spn = 1.0 # wing span
- lgt = 3.0 # channel half length
- hgt = 3.0 # channel half height
- wdt = 3.0 # channel width
- S_ref = 1.*spn # reference area
- c_ref = 1 # reference chord
- fms = 1.0 # farfield mesh size
- nms = 0.02 # nearfield mesh size
-
- # mesh the geometry
- print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
- tms['msh'].start()
- pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msLe' : nms, 'msTe' : nms, 'msF' : fms}
- msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
- tms['msh'].stop()
-
- # set the problem and add fluid, initial/boundary conditions
- tms['pre'].start()
- pbl, _, wk, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip', vsc=True)
- tms['pre'].stop()
-
- # solve problem
- config = {'nDim': dim, 'nSections':nSections, 'span':spn, 'nPoints':[100, 25], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 1e-8, 'degree': 1, 'neighbors': 50}
- iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw, config)
- vSolver = viscU.initBL(Re, M_inf, CFL0, nSections, 2)
- #iSolverAPI = viscAPI.dartAPI(floD.newton(pbl), bnd, wk, nSections, vInterp)
- coupler = viscC.Coupler(iSolverAPI, vSolver)
-
- coupler.Run()
-
- # display timers
- tms['total'].stop()
- print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
- print('CPU statistics')
- print(tms)
- print(coupler.tms)
-
- # visualize solution
- iSolverAPI.GetCp()
- iSolverAPI.iSolver.save(gmshWriter)
- floD.initViewer(pbl)
-
- # check results
- print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
- tests = CTests()
- if alpha == 3*math.pi/180 and M_inf == 0.3 and spn == 1.0:
- tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True))
- tests.add(CTest('CL', solver.Cl, 0.135, 5e-2))
- tests.add(CTest('CD', solver.Cd, 0.0062, 1e-2)) # Tranair (NF=0.0062, FF=0.0030), Panair 0.0035
- tests.add(CTest('CS', solver.Cs, 0.0121, 5e-2))
- tests.add(CTest('CM', solver.Cm, -0.0278, 1e-2))
- else:
- raise Exception('Test not defined for this flow')
- tests.run()
-
- # eof
- print('')
-
-if __name__ == "__main__":
- main()
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+## Compute lifting (linear or nonlinear) potential flow around a NACA 0012 rectangular wing
+# Adrien Crovato
+#
+# Test the nonlinear shock-capturing capability and the Kutta condition
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+import math
+import dart.utils as floU
+import dart.default as floD
+
+import vii.pyVII.vUtils as viscU
+import vii.pyVII.vCoupler2 as viscC
+import vii.pyVII.vInterpolator as vInterpol
+
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+def main():
+ # timer
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ # define flow variables
+ """Re = 11.72e6
+ alpha = 3.06*math.pi/180
+ M_inf = 0.8395"""
+ Re = 1e7
+ alpha = 0.*math.pi/180
+ M_inf = 0.
+ dim = 3
+ CFL0 = 1
+
+ nSections = 50
+
+ # define dimension and mesh size
+ spn = 1.0 # wing span
+ lgt = 3.0 # channel half length
+ hgt = 3.0 # channel half height
+ wdt = 3.0 # channel width
+ S_ref = 1.*spn # reference area
+ c_ref = 1 # reference chord
+ fms = 1.0 # farfield mesh size
+ nms = 0.02 # nearfield mesh size
+
+ # mesh the geometry
+ print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+ tms['msh'].start()
+ pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msLe' : nms, 'msTe' : nms, 'msF' : fms}
+ imsh = viscU.mesh('/home/paul/lab/dartflo/dart/models/n0012_3.geo', pars)
+ vmsh = viscU.mesh('/home/paul/lab/dartflo/vii/models/n0012_3_visc.geo', pars)
+ msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
+ tms['msh'].stop()
+
+ # set the problem and add fluid, initial/boundary conditions
+ tms['pre'].start()
+ # Replce tp = 'teTip' by tp = 'wakeTip' for oneraM6.geo
+ pbl, _, wk, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip', vsc=True)
+ tms['pre'].stop()
+
+ # solve problem
+ config = {'nDim': dim, 'nSections':nSections, 'span':spn, 'nPoints':[100, 25], 'eps':0.02, 'rbftype': 'linear', 'smoothing': 1e-6, 'degree': 0, 'neighbors': 50}
+ iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw, imsh, vmsh, config)
+ vSolver = viscU.initBL(Re, M_inf, CFL0, nSections, 2)
+ #iSolverAPI = viscAPI.dartAPI(floD.newton(pbl), bnd, wk, nSections, vInterp)
+ coupler = viscC.Coupler(iSolverAPI, vSolver)
+
+ coupler.Run()
+
+ # display timers
+ tms['total'].stop()
+ print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+ print('CPU statistics')
+ print(tms)
+ print(coupler.tms)
+
+ # visualize solution
+ iSolverAPI.GetCp()
+ iSolverAPI.iSolver.save(gmshWriter)
+ floD.initViewer(pbl)
+
+ # check results
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ if alpha == 3*math.pi/180 and M_inf == 0.3 and spn == 1.0:
+ tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True))
+ tests.add(CTest('CL', solver.Cl, 0.135, 5e-2))
+ tests.add(CTest('CD', solver.Cd, 0.0062, 1e-2)) # Tranair (NF=0.0062, FF=0.0030), Panair 0.0035
+ tests.add(CTest('CS', solver.Cs, 0.0121, 5e-2))
+ tests.add(CTest('CM', solver.Cm, -0.0278, 1e-2))
+ else:
+ raise Exception('Test not defined for this flow')
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
--
GitLab