From ebe605ce8e01ce28d51cdda4a289f44028acd4ba Mon Sep 17 00:00:00 2001
From: Paul <paul.dechamps@student.uliege.be>
Date: Mon, 11 Jul 2022 03:55:02 +0200
Subject: [PATCH] Half loop 3D
---
vii/models/n0012_3_visc.geo | 14 +--
vii/pyVII/vData.py | 166 +++++++++++++++++++++++++++++++++
vii/pyVII/vInterpolator.py | 177 +++++++++++++++++++++---------------
vii/tests/bli3.py | 6 +-
4 files changed, 281 insertions(+), 82 deletions(-)
create mode 100644 vii/pyVII/vData.py
diff --git a/vii/models/n0012_3_visc.geo b/vii/models/n0012_3_visc.geo
index b449ddc..5b4e93d 100644
--- a/vii/models/n0012_3_visc.geo
+++ b/vii/models/n0012_3_visc.geo
@@ -348,23 +348,23 @@ Surface(51) = {51};
//// MESHING ALGORITHM
//+
-nElmSpan = 50;
+nElmSpan = 40;
Transfinite Curve {65, 63, 62, 66, 61, 64, 241} = nElmSpan Using Progression 1;
-nElmChrodwise_wake = 20;
+nElmChrodwise_wake = 40;
//+
-Transfinite Curve {222, 221} = nElmChrodwise_wake Using Progression 1.01;
+Transfinite Curve {222, 221} = nElmChrodwise_wake Using Progression 1.05;
Transfinite Surface {51};
nChordwise_midWing = 40;
//+
-Transfinite Curve {8, 11, 5, 2, 122} = nChordwise_midWing Using Progression 1.01;
+Transfinite Curve {8, 11, 5, 2, 122} = nChordwise_midWing Using Progression 1;
//+
nChordwise_Le = 40;
-Transfinite Curve {9, 9, 10, 3, 4} = nChordwise_Le Using Progression 0.95;
+Transfinite Curve {-9, 10, -3, 4} = nChordwise_Le Using Progression 1.07;
//+
nChordwise_Te = 20;
-Transfinite Curve {7, 121, 12, 1, 6} = nChordwise_Te Using Progression 1.01;
+Transfinite Curve {7, 121, 12, 1, 6} = nChordwise_Te Using Progression 1;
Transfinite Surface {3};
Transfinite Surface {4};
@@ -394,4 +394,4 @@ Physical Surface("wing") = {1,2,3,11,12,13};
Physical Surface("wing_") = {4,5,6,14,15,16};
// Wake:
-Physical Surface("wake") = {51};
+Physical Surface("wake") = {51};
\ No newline at end of file
diff --git a/vii/pyVII/vData.py b/vii/pyVII/vData.py
new file mode 100644
index 0000000..0b44bf6
--- /dev/null
+++ b/vii/pyVII/vData.py
@@ -0,0 +1,166 @@
+#!/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.
+
+
+## VII Coupler
+# Paul Dechamps
+
+import fwk
+from fwk.coloring import ccolors
+
+from math import pi as M_PI
+import numpy as np
+
+class Data:
+ def __init__(self, msh, visc=0, PosSections= None):
+
+ self.Points, self.Cg = self.GetMesh(msh)
+
+ if visc:
+ self.Xp, self.Yp, self.Zp, self.sortedRows, self.Xc, self.Yc, self.Zc, self.sortedElems = self.SortMesh(self.Points, self.Cg)
+
+ # Creating sections
+ self.Sections = self.CreateSections(PosSections)
+
+ print(self.Sections)
+
+
+ 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:
+ alreadysaid.append(e.ptag.name)
+
+ 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))]
+ Xp = [np.zeros((nChord[i], 0)) for i in range(len(data))]
+ Yp = [np.zeros((nChord[i], 0)) for i in range(len(data))]
+ Zp = [np.zeros((nChord[i], 0)) for i in range(len(data))]
+
+ Xc = [np.zeros((nChord[i], 0)) for i in range(len(data))]
+ Yc = [np.zeros((nChord[i], 0)) for i in range(len(data))]
+ Zc = [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, :]
+ sort = sort[sort[:,0].argsort()]
+
+ sortedRows[iReg] = np.column_stack((sortedRows[iReg], sort[:,3]))
+ Xp[iReg] = np.column_stack((Xp[iReg], sort[:,0]))
+ Yp[iReg] = np.column_stack((Yp[iReg], sort[:,1]))
+ Zp[iReg] = np.column_stack((Zp[iReg], sort[:,2]))
+
+ 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, :]
+ sort = sort[sort[:,0].argsort()]
+
+ sortedElems[iReg] = np.column_stack((sortedElems[iReg], sort[:,3]))
+ Xc[iReg] = np.column_stack((Xc[iReg], sort[:,0]))
+ Yc[iReg] = np.column_stack((Yc[iReg], sort[:,1]))
+ Zc[iReg] = np.column_stack((Zc[iReg], sort[:,2]))
+
+ return Xp, Yp, Zp, sortedRows, Xc, Yc, Zc, sortedElems
+
+ def CreateSections(self, PosSections):
+ Sections = [[np.zeros((len(self.Xp[iReg]), 4)) for iReg in range(2)] for _ in range(len(PosSections))]
+
+ for iSec, y in enumerate(PosSections):
+ for iReg in range(2):
+ print(self.Yp[iReg][0,:])
+ print(np.where(self.Yp[iReg][0,:]==y)[0])
+ Sections[iSec][iReg][:,0] = self.Xp[iReg][:,np.where(self.Yp[iReg][0,:]==y)[0][0]]
+ Sections[iSec][iReg][:,1] = self.Yp[iReg][:,np.where(self.Yp[iReg][0,:]==y)[0][0]]
+ Sections[iSec][iReg][:,2] = self.Zp[iReg][:,np.where(self.Yp[iReg][0,:]==y)[0][0]]
+ Sections[iSec][iReg][:,3] = self.sortedRows[iReg][:,np.where(self.Yp[iReg][0,:]==y)[0][0]]
+
+ print(Sections[0][0])
+ return Sections
+
+
+
\ No newline at end of file
diff --git a/vii/pyVII/vInterpolator.py b/vii/pyVII/vInterpolator.py
index 573cc10..8f7210e 100644
--- a/vii/pyVII/vInterpolator.py
+++ b/vii/pyVII/vInterpolator.py
@@ -12,6 +12,8 @@ from scipy.interpolate import interp1d
from scipy.interpolate import RBFInterpolator
+from scipy.interpolate import interp2d
+
import fwk
@@ -25,7 +27,7 @@ class Interpolator:
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.Xp, self.Yp, self.Zp, self.sortedRows, self.Xc, self.Yc, self.Zc, self.sortedElems = self.SortMesh(self.vdata, self.vcg)
self.tms['Mesh sort'].stop()
print(self.tms)
@@ -61,22 +63,16 @@ class Interpolator:
# 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.stgPt=np.zeros(len(self.cfg['Sections']), dtype=int)
- self.xxPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(self.cfg['nSections'])]
+ self.xxPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(len(self.cfg['Sections']))]
#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'])]
+ self.DeltaStarPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(len(self.cfg['Sections']))]
+ self.Cg_Sec = [[np.zeros((len(self.Xc[iReg])-1, 3)) for iReg in range(2)] for _ in range(len(self.cfg['Sections']))]
def GetInvBC(self, vSolver, couplIter):
self.tms['InvToVisc'].start()
@@ -100,10 +96,6 @@ class Interpolator:
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])]
@@ -111,60 +103,88 @@ class Interpolator:
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()
+ for iSec, y in enumerate(self.cfg['Sections']):
- vSolver.SetVelocities(iSec, 0, UxUp, UyUp, UzUp)
- vSolver.SetVelocities(iSec, 1, UxLw, UyLw, UzLw)
- vSolver.SetVelocities(iSec, 2, Ux[1], Uy[1], Uz[1])
+ for iReg in range(2):
+
+ # Y Z inverted
+ x_pos = self.Xp[iReg][:,(abs(self.Yp[iReg][0,:]-y)).argmin()]
+ y_pos = self.Zp[iReg][:,(abs(self.Yp[iReg][0,:]-y)).argmin()]
+ z_pos = self.Yp[iReg][:,(abs(self.Yp[iReg][0,:]-y)).argmin()]
+ goodRows = self.sortedRows[iReg][:,(abs(self.Yp[iReg][0,:]-y)).argmin()]
+
+ ux = np.zeros(0)
+ uy = np.zeros(0)
+ uz = np.zeros(0)
+ me = np.zeros(0)
+ rhoe = np.zeros(0)
+ xxx = np.zeros(0)
+
+ for i in range(len(x_pos)):
+ # Y Z inverted
+ row = np.where(self.vdata[iReg][:,3]==goodRows[i])[0]
+ xxx = np.append(xxx, self.vdata[iReg][row,0])
+ ux = np.append(ux, Ux[iReg][row])
+ uy = np.append(uy, Uz[iReg][row])
+ uz = np.append(uz, Uy[iReg][row])
+ me = np.append(me, Me[iReg][row])
+ rhoe = np.append(rhoe, Rhoe[iReg][row])
+
+ if couplIter == 0:
+
+ if iReg == 0:
+ self.stgPt[iSec] = np.argmin(ux**2+uy**2)
+ xUp, xLw = self.__Remesh(abs(x_pos), self.stgPt[iSec])
+ yUp, yLw = self.__Remesh(y_pos, self.stgPt[iSec])
+ zUp, zLw = self.__Remesh(z_pos, self.stgPt[iSec])
+
+ vSolver.SetGlobMesh(iSec, 0, xUp, yUp, zUp)
+ vSolver.SetGlobMesh(iSec, 1, xLw, yLw, zLw)
+
+ self.DeltaStarPrev[iSec][0] = np.zeros(len(xUp))
+ self.DeltaStarPrev[iSec][1] = np.zeros(len(xLw))
+
+ self.xxPrev[iSec][0] = np.zeros(len(xUp))
+ self.xxPrev[iSec][1] = np.zeros(len(xLw))
+ else:
+ vSolver.SetGlobMesh(iSec, 2, x_pos, y_pos, z_pos)
+ self.DeltaStarPrev[iSec][2] = np.zeros(len(x_pos))
+ self.xxPrev[iSec][2] = np.zeros(len(x_pos))
+
+ for i in range(len(x_pos)-1):
+ self.Cg_Sec[iSec][iReg][i,0] = x_pos[i+1] + x_pos[i]
+ self.Cg_Sec[iSec][iReg][i,1] = y_pos[i+1] + y_pos[i]
+ self.Cg_Sec[iSec][iReg][i,2] = z_pos[i+1] + z_pos[i]
+ self.Cg_Sec[iSec][iReg]/=2
+
+ plt.plot(x_pos, y_pos, 'x-')
+ plt.plot(self.Cg_Sec[iSec][iReg][:,0], self.Cg_Sec[iSec][iReg][:,1], 'o')
+ plt.show()
+
+ if iReg == 0:
+ UxUp, UxLw = self.__Remesh(ux, self.stgPt[iSec])
+ UyUp, UyLw = self.__Remesh(uy, self.stgPt[iSec])
+ UzUp, UzLw = self.__Remesh(uz, self.stgPt[iSec])
- MeUp, MeLw = self.__Remesh(Me[0], self.stgPt[iSec])
+ vSolver.SetVelocities(iSec, 0, UxUp, UyUp, UzUp)
+ vSolver.SetVelocities(iSec, 1, UxLw, UyLw, UzLw)
- vSolver.SetMe(iSec, 0, MeUp)
- vSolver.SetMe(iSec, 1, MeLw)
- vSolver.SetMe(iSec, 2, Me[1])
+ MeUp, MeLw = self.__Remesh(me, self.stgPt[iSec])
+ vSolver.SetMe(iSec, 0, MeUp)
+ vSolver.SetMe(iSec, 1, MeLw)
- RhoUp, RhoLw = self.__Remesh(Rhoe[0], self.stgPt[iSec])
+ RhoUp, RhoLw = self.__Remesh(rhoe, self.stgPt[iSec])
+ vSolver.SetRhoe(iSec, 0, RhoUp)
+ vSolver.SetRhoe(iSec, 1, RhoLw)
+ else:
- vSolver.SetRhoe(iSec, 0, RhoUp)
- vSolver.SetRhoe(iSec, 1, RhoLw)
- vSolver.SetRhoe(iSec, 2, Rhoe[1])
+ vSolver.SetVelocities(iSec, 2, ux, uy, uz)
+ vSolver.SetMe(iSec, 2, me)
+ vSolver.SetRhoe(iSec, 2, rhoe)
- for iReg in range(3):
- vSolver.SetDeltaStarExt(iSec, iReg, self.DeltaStarPrev[iSec][iReg])
- vSolver.SetxxExt(iSec, iReg, self.xxPrev[iSec][iReg])
-
-
+ 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):
@@ -173,20 +193,19 @@ class Interpolator:
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)))
+ blwAll = interp2d(self.Cg_Sec[0])
+
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(x[0][:,0], x[0][:,2], blw[0], color = 'red')
@@ -590,9 +609,7 @@ class Interpolator:
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
@@ -619,6 +636,9 @@ class Interpolator:
cg_pt[0] = -cg_pt[0]
cg[0] = np.vstack((cg[0], cg_pt))
+ print(msh.elems.size())
+ print(np.shape(cg[0]))
+ print(np.shape(data[0]))
return data, cg
def SortMesh(self, data, cg):
@@ -639,6 +659,9 @@ class Interpolator:
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))]
+ Xp = [np.zeros((nChord[i], 0)) for i in range(len(data))]
+ Yp = [np.zeros((nChord[i], 0)) for i in range(len(data))]
+ Zp = [np.zeros((nChord[i], 0)) for i in range(len(data))]
for iReg, reg in enumerate(data):
@@ -650,8 +673,12 @@ class Interpolator:
sections[iy] = reg[idx,1]
sort = reg[abs(reg[:,1]-sections[iy])<1e-3, :]
+ sort = sort[sort[:,0].argsort()]
sortedRows[iReg] = np.column_stack((sortedRows[iReg], sort[:,3]))
+ Xp[iReg] = np.column_stack((Xp[iReg], sort[:,0]))
+ Yp[iReg] = np.column_stack((Yp[iReg], sort[:,1]))
+ Zp[iReg] = np.column_stack((Zp[iReg], sort[:,2]))
sections_cg = cg[0][:,1]
for i in range(len(sections_cg)):
@@ -660,6 +687,10 @@ class Interpolator:
nChord_cg = [len(cg[i][abs(cg[i][:,1]-sections[0])<1e-3, 0]) for i in range(len(cg))]
+ Xc = [np.zeros((nChord_cg[i], 0)) for i in range(len(data))]
+ Yc = [np.zeros((nChord_cg[i], 0)) for i in range(len(data))]
+ Zc = [np.zeros((nChord_cg[i], 0)) for i in range(len(data))]
+
sortedElems = [np.zeros((nChord_cg[i], 0)) for i in range(len(cg))]
for iReg, reg in enumerate(cg):
@@ -672,10 +703,14 @@ class Interpolator:
sections_cg[iy] = reg[idx,1]
sort = reg[abs(reg[:,1]-sections_cg[iy])<1e-3, :]
-
+ sort = sort[sort[:,0].argsort()]
+
sortedElems[iReg] = np.column_stack((sortedElems[iReg], sort[:,3]))
+ Xc[iReg] = np.column_stack((Xc[iReg], sort[:,0]))
+ Yc[iReg] = np.column_stack((Yc[iReg], sort[:,1]))
+ Zc[iReg] = np.column_stack((Zc[iReg], sort[:,2]))
- return sortedRows, sortedElems
+ return Xp, Yp, Zp, sortedRows, Xc, Yc, Zc, sortedElems
diff --git a/vii/tests/bli3.py b/vii/tests/bli3.py
index 29c99f8..139f34e 100644
--- a/vii/tests/bli3.py
+++ b/vii/tests/bli3.py
@@ -52,8 +52,6 @@ def main():
dim = 3
CFL0 = 1
- nSections = 50
-
# define dimension and mesh size
spn = 1.0 # wing span
lgt = 3.0 # channel half length
@@ -80,9 +78,9 @@ def main():
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}
+ config = {'nDim': dim, 'Sections':[0.02, 0.041], 'span':spn, 'rbftype': 'linear', 'smoothing': 1e-8, 'degree': 0, 'neighbors': 5}
iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw, imsh, vmsh, config)
- vSolver = viscU.initBL(Re, M_inf, CFL0, nSections, 2)
+ vSolver = viscU.initBL(Re, M_inf, CFL0, len(config['Sections']), 2)
#iSolverAPI = viscAPI.dartAPI(floD.newton(pbl), bnd, wk, nSections, vInterp)
coupler = viscC.Coupler(iSolverAPI, vSolver)
--
GitLab