diff --git a/dart/pyVII/vAPI3.py b/dart/pyVII/vAPI3.py
new file mode 100644
index 0000000000000000000000000000000000000000..b5718a536de65cc63a8cd36364b851f344969108
--- /dev/null
+++ b/dart/pyVII/vAPI3.py
@@ -0,0 +1,107 @@
+import numpy as np
+import dart.pyVII.vStruct as vStruct
+from matplotlib import pyplot as plt
+from scipy import interpolate
+
+class dartAPI:
+
+ def __init__(self, _dartSolver, bnd_wing, bnd_wake, _nSections, _interp):
+ self.nSections = _nSections
+ self.iSolver = _dartSolver
+ self.bndWing = bnd_wing
+ self.bndWk = bnd_wake
+ #self.invStruct = [[vStruct.Airfoil(bnd_airfoil), vStruct.Wake(bnd_wake)] for _ in range(self.nSections)] # Parameters passing protocols
+ self.viscStruct = [[vStruct.viscousIPS("upper"), vStruct.viscousIPS("lower"), vStruct.viscousIPS("wake")] for _ in range(self.nSections)]
+ self.interp = _interp
+
+ def GetInvBC(self, vSolver, couplIter):
+ Ux = np.zeros(len(self.bndWing.nodes))
+ Uy = np.zeros(len(self.bndWing.nodes))
+ Uz = np.zeros(len(self.bndWing.nodes))
+ Me = np.zeros(len(self.bndWing.nodes))
+ Rho = np.zeros(len(self.bndWing.nodes))
+ UxWk = np.zeros(len(self.bndWk.nodMap))
+ UyWk = np.zeros(len(self.bndWk.nodMap))
+ UzWk = np.zeros(len(self.bndWk.nodMap))
+ MeWk = np.zeros(len(self.bndWk.nodMap))
+ RhoWk = np.zeros(len(self.bndWk.nodMap))
+
+ # From wing
+ for node in range(len(self.bndWing.nodes)):
+ Ux[node] = self.iSolver.U[self.bndWing.nodes[node].row][0]
+ Uy[node] = self.iSolver.U[self.bndWing.nodes[node].row][1]
+ Uz[node] = self.iSolver.U[self.bndWing.nodes[node].row][2]
+
+ Me[node] = self.iSolver.M[self.bndWing.nodes[node].row]
+ Rho[node] = self.iSolver.rho[self.bndWing.nodes[node].row]
+
+ # From wake
+ for iNodeWk, n in enumerate(self.bndWk.nodMap):
+ UxWk[iNodeWk] = self.iSolver.U[n.row][0]
+ UyWk[iNodeWk] = self.iSolver.U[n.row][1]
+ UzWk[iNodeWk] = self.iSolver.U[n.row][2]
+
+ MeWk[iNodeWk] = self.iSolver.M[n.row]
+ RhoWk[iNodeWk] = self.iSolver.rho[n.row]
+
+ mapedUx = self.interp.InterpToSections(Ux, UxWk)
+ mapedUy = self.interp.InterpToSections(Uy, UyWk)
+ mapedUz = self.interp.InterpToSections(Uz, UzWk)
+
+ """plt.plot(self.interp.Sections[3][0][:,0], np.sqrt(mapedUx[3][0]**2 + mapedUy[3][0]**2 +mapedUz[3][0]**2))
+ plt.plot(self.interp.Sections[3][1][:,0], np.sqrt(mapedUx[3][1]**2 + mapedUy[3][1]**2 +mapedUz[3][1]**2))
+ plt.plot(self.interp.Sections[3][2][:,0], np.sqrt(mapedUx[3][2]**2 + mapedUy[3][2]**2 +mapedUz[3][2]**2))
+ plt.show()"""
+
+ mapedMe = self.interp.InterpToSections(Me, MeWk)
+ mapedRho = self.interp.InterpToSections(Rho, RhoWk)
+
+ for iSec in range(len(self.interp.Sections)):
+ for iReg in range(3):
+ vSolver.SetGlobMesh(iSec, iReg, self.interp.Sections[iSec][iReg][:,0],
+ self.interp.Sections[iSec][iReg][:,1],
+ self.interp.Sections[iSec][iReg][:,2])
+ vSolver.SetVelocities(iSec, iReg, mapedUx[iSec][iReg], mapedUy[iSec][iReg], mapedUz[iSec][iReg])
+ vSolver.SetMe(iSec, iReg, mapedMe[iSec][iReg])
+ vSolver.SetRhoe(iSec, iReg, mapedRho[iSec][iReg])
+ if couplIter == 0:
+ self.viscStruct[iSec][iReg].xxExt = np.zeros(len(self.interp.Sections[iSec][iReg][:,0]))
+ self.viscStruct[iSec][iReg].DeltaStarExt = np.zeros(len(self.interp.Sections[iSec][iReg][:,0]))
+ vSolver.SetxxExt(iSec, iReg, self.viscStruct[iSec][iReg].xxExt)
+ vSolver.SetDeltaStarExt(iSec, iReg, self.viscStruct[iSec][iReg].DeltaStarExt)
+
+ def SetBlowingVel(self, vSolver):
+ """ Collects blowing velocities from inviscid solver, maps them to inviscid mesh and """
+
+ for iSec in range(self.nSections):
+ for iRegion in range(len(self.viscStruct[iSec])):
+ self.viscStruct[iSec][iRegion].BlowingVel = np.asarray(vSolver.ExtractBlowingVel(iSec, iRegion))
+ self.viscStruct[iSec][iRegion].DeltaStarExt = vSolver.ExtractDeltaStar(iSec, iRegion)
+ self.viscStruct[iSec][iRegion].xxExt = vSolver.Extractxx(iSec, iRegion)
+
+ self.interp.InterpToElements(self.viscStruct)
+
+ 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
+
\ No newline at end of file
diff --git a/dart/pyVII/vCoupler.py b/dart/pyVII/vCoupler.py
index 2b69eca0b9b50969d57677bf2eba9190f4b2615e..8c623cab11c7408867949ee15df5b479450b1bb4 100644
--- a/dart/pyVII/vCoupler.py
+++ b/dart/pyVII/vCoupler.py
@@ -27,6 +27,7 @@ import math
import fwk
import dart.pyVII.vUtils as viscU
import numpy as np
+from matplotlib import pyplot as plt
class Coupler:
def __init__(self, _iSolverAPI, vSolver) -> None:
@@ -77,7 +78,7 @@ class Coupler:
error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
- if error <= self.couplTol and exitCode==0:
+ if error <= self.couplTol:
print(ccolors.ANSI_GREEN,'{:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}\n'.format(couplIter,
self.iSolverAPI.GetCl(), self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)),ccolors.ANSI_RESET)
@@ -119,82 +120,4 @@ class Coupler:
coeffTable[1,ind] = self.iSolver.Cl
coeffTable[2,ind] = self.vSolver.Cdt
- return coeffTable
-
- def __ImposeBlwVel(self):
-
- # Retreive and set blowing velocities.
-
- for iRegion in range(3):
- self.blwVel[iRegion] = np.asarray(self.vSolver.ExtractBlowingVel(iRegion))
- self.deltaStar[iRegion] = np.asarray(self.vSolver.ExtractDeltaStar(iRegion))
- self.xx[iRegion] = np.asarray(self.vSolver.ExtractXx(iRegion))
-
- blwVelAF = np.concatenate((self.blwVel[0], self.blwVel[1]))
- deltaStarAF = np.concatenate((self.deltaStar[0], self.deltaStar[1][1:])) # Remove stagnation point doublon.
- xxAF = np.concatenate((self.xx[0], self.xx[1][1:])) # Remove stagnation point doublon.
- blwVelWK = self.blwVel[2]
- deltaStarWK = self.deltaStar[2]
- xxWK = self.xx[2]
-
- self.group[0].u = blwVelAF[self.group[0].connectListElems.argsort()]
- self.group[1].u = blwVelWK[self.group[1].connectListElems.argsort()]
-
- self.group[0].deltaStar = deltaStarAF[self.group[0].connectListNodes.argsort()]
- self.group[1].deltaStar = deltaStarWK[self.group[1].connectListNodes.argsort()]
-
- self.group[0].xx = xxAF[self.group[0].connectListNodes.argsort()]
- self.group[1].xx = xxWK[self.group[1].connectListNodes.argsort()]
-
- for n in range(0, len(self.group)):
- for i in range(0, self.group[n].nE):
- self.group[n].boundary.setU(i, self.group[n].u[i])
-
- def __ImposeInvBC(self, couplIter):
-
- for n in range(0, len(self.group)):
-
- for i in range(0, len(self.group[n].boundary.nodes)):
-
- self.group[n].v[i,0] = self.iSolver.U[self.group[n].boundary.nodes[i].row][0]
- self.group[n].v[i,1] = self.iSolver.U[self.group[n].boundary.nodes[i].row][1]
- self.group[n].v[i,2] = self.iSolver.U[self.group[n].boundary.nodes[i].row][2]
- self.group[n].Me[i] = self.iSolver.M[self.group[n].boundary.nodes[i].row]
- self.group[n].rhoe[i] = self.iSolver.rho[self.group[n].boundary.nodes[i].row]
-
- dataIn=[None,None]
- for n in range(0,len(self.group)):
- self.group[n].connectListNodes, self.group[n].connectListElems, dataIn[n] = self.group[n].connectList()
- fMeshDict, _, _ = GroupConstructor().mergeVectors(dataIn, 0, 1)
-
- if ((len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]]) != self.nElmUpperPrev) or couplIter==0):
-
- # print('STAG POINT MOVED')
-
- # Initialize mesh
- self.vSolver.InitMeshes(0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]], fMeshDict['globCoord'][:fMeshDict['bNodes'][1]])
- self.vSolver.InitMeshes(1, fMeshDict['locCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
- self.vSolver.InitMeshes(2, fMeshDict['locCoord'][fMeshDict['bNodes'][2]:], fMeshDict['globCoord'][fMeshDict['bNodes'][2]:])
-
- self.vSolver.SendExtVariables(0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]], fMeshDict['deltaStarExt'][:fMeshDict['bNodes'][1]])
- self.vSolver.SendExtVariables(1, fMeshDict['locCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
- self.vSolver.SendExtVariables(2, fMeshDict['locCoord'][fMeshDict['bNodes'][2]:], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][2]:])
-
- else:
-
- self.vSolver.SendExtVariables(0, fMeshDict['xxExt'][:fMeshDict['bNodes'][1]], fMeshDict['deltaStarExt'][:fMeshDict['bNodes'][1]])
- self.vSolver.SendExtVariables(1, fMeshDict['xxExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
- self.vSolver.SendExtVariables(2, fMeshDict['xxExt'][fMeshDict['bNodes'][2]:], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][2]:])
-
- self.nElmUpperPrev = len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]])
- self.vSolver.SendInvBC(0, fMeshDict['vtExt'][:fMeshDict['bNodes'][1]], fMeshDict['Me'][:fMeshDict['bNodes'][1]], fMeshDict['rho'][:fMeshDict['bNodes'][1]])
- self.vSolver.SendInvBC(1, fMeshDict['vtExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['Me'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['rho'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
- self.vSolver.SendInvBC(2, fMeshDict['vtExt'][fMeshDict['bNodes'][2]:], fMeshDict['Me'][fMeshDict['bNodes'][2]:], fMeshDict['rho'][fMeshDict['bNodes'][2]:])
- # print('Inviscid boundary imposed')
-
-
-
-
-
-
-
+ return coeffTable
\ No newline at end of file
diff --git a/dart/pyVII/vInterpolator3.py b/dart/pyVII/vInterpolator3.py
new file mode 100644
index 0000000000000000000000000000000000000000..c02a5384c6ee64bcde555d5d6725c62325021ef2
--- /dev/null
+++ b/dart/pyVII/vInterpolator3.py
@@ -0,0 +1,195 @@
+import numpy as np
+import math as m
+from matplotlib import pyplot as plt
+from scipy.interpolate import bisplrep
+from scipy.interpolate import bisplev
+from scipy import stats
+from scipy.interpolate import RBFInterpolator
+
+import sys
+
+class Interpolator:
+
+ def __init__(self, bnd, wk, nSections, blw) -> None:
+ self.bndWing = bnd
+ self.bndWake = wk
+ self.blwWing = blw[0]
+ self.blwWake = blw[1]
+ self.Sections, self.wing, self.wake = self.GetWing(bnd, wk, nSections)
+
+
+ def GetWing(self, bnd, wk, nSections, nPoints=[500, 25], eps=5e-2, k=[5, 5]):
+
+ "Extract and sort the node coordinates of the airfoil."
+ wing = np.zeros((len(bnd.nodes),3))
+ wake = np.zeros((len(wk.nodMap),3))
+
+ # Extract coordinates.
+ for iNode, n in enumerate(bnd.nodes):
+ for iDir in range(len(wing[0,:])):
+ wing[iNode,iDir] = n.pos[iDir]
+
+ for iNode, n in enumerate(wk.nodMap):
+ for iDir in range(len(wake[0,:])):
+ wake[iNode, iDir] = n.pos[iDir]
+
+ # Separate upper and lower sides
+
+ wingUp = wing[wing[:,2]>=0]
+ wingLw = wing[wing[:,2]<=0]
+
+ wingUp = np.delete(wingUp, wingUp[:,1]>1, 0)
+ wingLw = np.delete(wingLw, wingLw[:,1]>1, 0)
+
+ xUp = wingUp[:,0].copy()
+ xLw = wingLw[:,0].copy()
+
+ np.set_printoptions(threshold=sys.maxsize)
+
+ sUp = (len(xUp)-np.sqrt(2*len(xUp)) + len(xUp)+np.sqrt(2*len(xUp)))/2
+ sLw = (len(xLw)-np.sqrt(2*len(xLw)) + len(xLw)+np.sqrt(2*len(xLw)))/2
+
+ tckUp = bisplrep(wingUp[:,0], wingUp[:,1], wingUp[:,2], kx=k[0], ky=k[1], s=sUp)
+ tckLw = bisplrep(wingLw[:,0], wingLw[:,1], wingLw[:,2], kx=k[0], ky=k[1], s=sLw)
+
+ yDistri = np.linspace(min(wing[:,1]), max(wing[:,1]), nSections)
+
+ Sections = [[np.zeros((nPoints[1] if i==2 else nPoints[0], 3)) for i in range(3)] for _ in range(nSections)]
+
+ for iSec, y in enumerate(yDistri):
+ for iReg in range(3):
+ if iReg == 0:
+ Sections[iSec][iReg][:,0] = np.linspace(min(wing[abs(wing[:,1]-y)<eps, 0]), max(wing[abs(wing[:,1]-y)<eps, 0]), nPoints[0])
+ elif iReg == 1:
+ Sections[iSec][iReg][:,0] = np.linspace(min(wing[abs(wing[:,1]-y)<eps, 0]), max(wing[abs(wing[:,1]-y)<eps, 0]), nPoints[0])
+ elif iReg == 2:
+ print('sec', iSec, 'min', min(wake[:, 0]), 'max', max(wake[:, 0]))
+ Sections[iSec][iReg][:,0] = np.linspace(max(wing[abs(wing[:,1]-y)<eps, 0]), max(wake[:, 0]), nPoints[1])
+ Sections[iSec][iReg][:,1] = y
+ for iPt in range(len(Sections[iSec][iReg][:,0])):
+ if iReg==0:
+ Sections[iSec][iReg][iPt,2] = bisplev(Sections[iSec][iReg][iPt, 0], Sections[iSec][iReg][iPt, 1], tckUp)
+ elif iReg==1:
+ Sections[iSec][iReg][iPt,2] = bisplev(Sections[iSec][iReg][iPt, 0], Sections[iSec][iReg][iPt, 1], tckLw)
+ elif iReg == 2:
+ Sections[iSec][iReg][iPt,2] = 0
+ else:
+ raise RuntimeError("dart::vInterpolator3 l.83 0<=iReg<=2")
+
+ # Add leading edge point
+
+ for iSec in range(len(Sections)):
+ for iReg in range(2):
+ Sections[iSec][iReg][0,2] = 0
+
+ """fig = plt.figure()
+ ax = fig.add_subplot(projection='3d')
+ ax.scatter(wing[:,0], wing[:,1], wing[:,2], s=0.3)
+ ax.scatter(wake[:,0], wake[:,1], wake[:,2], s=0.3)
+ for iSec in range(len(Sections)):
+ ax.scatter(Sections[iSec][0][:,0], Sections[iSec][0][:,1], Sections[iSec][0][:,2], color='red')
+ ax.scatter(Sections[iSec][1][:,0], Sections[iSec][1][:,1], Sections[iSec][1][:,2], color='red')
+ ax.scatter(Sections[iSec][2][:,0], Sections[iSec][2][:,1], Sections[iSec][2][:,2], color='red')
+ ax.set_zlim([-0.5, 0.5])
+ ax.set_xlabel('x')
+ ax.set_ylabel('y')
+ ax.set_zlabel('z')
+ plt.show()
+
+ plt.plot(Sections[0][0][:,0], Sections[0][0][:,2], 'x-', color='red')
+ plt.plot(Sections[0][1][:,0], Sections[0][1][:,2], 'x-', color='red')
+ plt.plot(wingUp[:,0], wingUp[:,2], 'o', color='blue')
+ plt.plot(wingLw[:,0], wingLw[:,2], 'o', color='blue')
+ plt.ylim([-0.5, 0.5])
+ plt.show()"""
+ return Sections, wing, wake
+
+ def InterpToSections(self, var, varWk):
+
+ varWk = np.reshape(varWk, [len(varWk), 1])
+ completeWake = np.hstack((self.wake, varWk))
+ _, idx = np.unique(completeWake[:,:3], axis=0, return_index=True)
+ completeWake = completeWake[np.sort(idx)]
+
+ var = np.reshape(var, [len(var), 1])
+ completeWing = np.hstack((self.wing, var))
+
+ wingUp = completeWing[completeWing[:,2]>=0]
+ wingLw = completeWing[completeWing[:,2]<=0]
+
+ mappedVar = [[np.zeros(0) for _ in range(3)] for _ in range(len(self.Sections))]
+
+ for iSec in range(len(mappedVar)):
+ for iReg in range(len(mappedVar[iSec])):
+ if iReg==0:
+ mappedVar[iSec][iReg] = RBFInterpolator(wingUp[:,:3], wingUp[:,3], neighbors=75, smoothing=1.5)(self.Sections[iSec][iReg])
+ if iReg==1:
+ mappedVar[iSec][iReg] = RBFInterpolator(wingLw[:,:3], wingLw[:,3], neighbors=75, smoothing=1.5)(self.Sections[iSec][iReg])
+ if iReg==2:
+ mappedVar[iSec][iReg] = RBFInterpolator(completeWake[:,:2], completeWake[:,3], smoothing=1.5)(self.Sections[iSec][iReg][:,:2])
+ return mappedVar
+
+ def InterpToElements(self, viscStruct):
+
+ blwwing = np.zeros((len(self.blwWing.nodes),3))
+ blwwake = np.zeros((len(self.blwWake.nodes),3))
+
+ # Extract coordinates.
+ for iNode, n in enumerate(self.blwWing.nodes):
+ for iDir in range(len(blwwing[0,:])):
+ blwwing[iNode,iDir] = n.pos[iDir]
+
+ for iNode, n in enumerate(self.blwWake.nodes):
+ for iDir in range(len(blwwake[0,:])):
+ blwwake[iNode, iDir] = n.pos[iDir]
+
+ coord_Up=np.zeros((3))
+ coord_Lw=np.zeros((3))
+ coord_Wk=np.zeros((3))
+ BlwTot_Up=np.zeros(1)
+ BlwTot_Lw=np.zeros(1)
+ BlwTot_Wk=np.zeros(1)
+ for iSec in range(len(viscStruct)):
+
+ coord_Up = np.vstack((coord_Up, self.Sections[iSec][0]))
+ if iSec!=0:
+ BlwTot_Up = np.insert(BlwTot_Up, -1, 0)
+ BlwTot_Up = np.concatenate([BlwTot_Up, viscStruct[iSec][0].BlowingVel], axis=None)
+
+ coord_Lw = np.vstack((coord_Lw, self.Sections[iSec][1]))
+ if iSec!=0:
+ BlwTot_Lw = np.insert(BlwTot_Lw, -1, 0)
+ BlwTot_Lw = np.concatenate([BlwTot_Lw, viscStruct[iSec][1].BlowingVel], axis=None)
+
+ coord_Wk = np.vstack((coord_Wk, self.Sections[iSec][2]))
+ if iSec!=0:
+ BlwTot_Wk = np.insert(BlwTot_Wk, -1, 0)
+ BlwTot_Wk = np.concatenate([BlwTot_Wk, viscStruct[iSec][2].BlowingVel], axis=None)
+
+ coord_Up = np.delete(coord_Up, 0, axis=0)
+ coord_Lw = np.delete(coord_Lw, 0, axis=0)
+ coord_Wk = np.delete(coord_Wk, 0, axis=0)
+
+ coord = np.vstack((coord_Up, coord_Lw))
+ blw = np.concatenate((BlwTot_Up, BlwTot_Lw), axis=0)
+
+
+ blw = np.reshape(blw, [len(blw), 1])
+ BlwTot_Wk = np.reshape(BlwTot_Wk, [len(BlwTot_Wk), 1])
+ completeWing = np.hstack((coord, blw))
+ completeWake = np.hstack((coord_Wk, BlwTot_Wk))
+
+ _, idx = np.unique(completeWing[:,:3], axis=0, return_index=True)
+ completeWing = completeWing[np.sort(idx)]
+
+ _, idxWk = np.unique(completeWake[:,:3], axis=0, return_index=True)
+ completeWake = completeWake[np.sort(idxWk)]
+
+ mappedBlowingWing = RBFInterpolator(completeWing[:,:3], completeWing[:,3], smoothing=1.5)(blwwing)
+ mappedBlowingWake = RBFInterpolator(completeWake[:,:2], completeWake[:,3], smoothing=1.5)(blwwake[:,:2])
+
+
+ for iNode in range(len(blwwing)):
+ self.blwWing.setU(iNode, mappedBlowingWing[iNode])
+ for iNode in range(len(blwwake)):
+ self.blwWake.setU(iNode, mappedBlowingWake[iNode])
\ No newline at end of file
diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index 60b71287e7ffc320e5e44832197094d9cb826dd9..c1ed9ac5717d3cb311712ce3fb061a8eb8917baf 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -20,7 +20,7 @@ BLRegion::~BLRegion()
{
delete closSolver;
delete mesh;
- std::cout << "~BLRegion()\n";
+ // std::cout << "~BLRegion()\n";
} // end ~BLRegion
void BLRegion::SetMesh(std::vector<double> x,
diff --git a/dart/src/wClosures.cpp b/dart/src/wClosures.cpp
index edaf5109f34b53add692ab6cdd99bd8c563ce5ef..4a91354875a4eee6dca9804a71d83bf94e5c7c03 100644
--- a/dart/src/wClosures.cpp
+++ b/dart/src/wClosures.cpp
@@ -11,7 +11,7 @@ Closures::Closures() {}
Closures::~Closures()
{
- std::cout << "~Closures()\n";
+ // std::cout << "~Closures()\n";
} // end ~Closures
std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret, double Me, std::string name)
diff --git a/dart/src/wProblem.cpp b/dart/src/wProblem.cpp
index 73d6f912234591de19674f8ef6cd03e9275d95ee..164dd47af954f2690e156752fafb42e106675406 100644
--- a/dart/src/wProblem.cpp
+++ b/dart/src/wProblem.cpp
@@ -253,8 +253,8 @@ void Problem::check() const
}
}
// Blowing B.C.
- if (!bBCs.empty())
- throw std::runtime_error("Blowing boundary conditions are not supported for 3D problems!\n");
+ /*if (!bBCs.empty())
+ throw std::runtime_error("Blowing boundary conditions are not supported for 3D problems!\n");*/
}
// Two-dimension problem
else if (nDim == 2)
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index e2d6e91a381a0dc6e59bbe00b44193f9dbfd04b6..3e993035f3b4b33fa7fd1512e04bdf93211ed2bd 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -95,12 +95,14 @@ ViscSolver::~ViscSolver()
*/
int ViscSolver::Run(unsigned int couplIter)
{
+ std::cout << "Starting solve" << std::endl;
bool lockTrans;
int solverExitCode = 0;
int pointExitCode;
for (size_t iSec=0; iSec<Sections.size(); ++iSec)
{
+ std::cout << "Section: " << iSec << std::endl;
/* Prepare time integration */
@@ -272,7 +274,6 @@ int ViscSolver::Run(unsigned int couplIter)
} // End for iPoints
Sections[iSec][iRegion]->ComputeBlwVel();
} // End for iRegion
-
} // End for iSections
ComputeDrag(Sections[0]);
@@ -433,7 +434,6 @@ void ViscSolver::ComputeDrag(std::vector<BLRegion *> bl)
void ViscSolver::SetGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z)
{
- std::cout << "len x " << _x.size() << std::endl;
Sections[iSec][iRegion]->SetMesh(_x, _y, _z);
}
@@ -446,7 +446,9 @@ void ViscSolver::SetVelocities(size_t iSec, size_t iRegion, std::vector<double>
{
alpha = std::atan2(Sections[iSec][iRegion]->mesh->Gety(iPoint+1) - Sections[iSec][iRegion]->mesh->Gety(iPoint),
Sections[iSec][iRegion]->mesh->Getx(iPoint+1) - Sections[iSec][iRegion]->mesh->Getx(iPoint));
+
Vt[iPoint] = vx[iPoint]*std::cos(alpha) + vy[iPoint]*std::sin(alpha);
+ // std::cout << iSec << iRegion << iPoint << "vx " << vx[iPoint] << "vy" << vy[iPoint] << "vz" << vz[iPoint] << "vt" << Vt[iPoint] << std::endl;
}
Vt.back() = vx.back()*cos(alpha) + vy.back()*std::sin(alpha);
@@ -465,13 +467,11 @@ void ViscSolver::SetRhoe(size_t iSec, size_t iRegion, std::vector<double> _Rhoe)
void ViscSolver::SetxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt)
{
- std::cout << "len xxExt " << _xxExt.size() << std::endl;
Sections[iSec][iRegion]->mesh->SetExt(_xxExt);
}
void ViscSolver::SetDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double> _DeltaStarExt)
{
- std::cout << "len DeltaStarExt " << _DeltaStarExt.size() << std::endl;
Sections[iSec][iRegion]->invBnd->SetDeltaStar(_DeltaStarExt);
}
diff --git a/dart/tests/bli3.py b/dart/tests/bli3.py
new file mode 100644
index 0000000000000000000000000000000000000000..17e63da7ee8717b753b528fe19310260d0149bd0
--- /dev/null
+++ b/dart/tests/bli3.py
@@ -0,0 +1,111 @@
+#!/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.vAPI3 as viscAPI
+import dart.pyVII.vCoupler as viscC
+import dart.pyVII.vInterpolator3 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 = 1e7
+ alpha = 3*math.pi/180
+ M_inf = 0.3
+ 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
+ vInterp = vInterpol.Interpolator(bnd, wk, nSections, blw)
+ vSolver = viscU.initBL(Re, M_inf, CFL0, nSections)
+ 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)
+
+ # visualize solution
+ 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()
diff --git a/dart/tests/bli_lowLift.py b/dart/tests/bli_lowLift.py
index 1494838939216721161c747968ab7c32bcfdc4aa..3c79f984844fbba1040784e1c55a1c13417651b8 100644
--- a/dart/tests/bli_lowLift.py
+++ b/dart/tests/bli_lowLift.py
@@ -97,7 +97,6 @@ def main():
tms['solver'].start()
vInterp = vInterpol.Interpolator(bnd)
- quit()
vSolver = viscU.initBL(Re, M_inf, CFL0, nSections)
iSolverAPI = viscAPI.dartAPI(floD.newton(pbl), blw[0], blw[1], nSections)
coupler = viscC.Coupler(iSolverAPI, vSolver)