diff --git a/blast/interfaces/blData.py b/blast/interfaces/blData.py new file mode 100644 index 0000000000000000000000000000000000000000..d2779f4bdb3a6e35d6422fb65f1d511dbbdda184 --- /dev/null +++ b/blast/interfaces/blData.py @@ -0,0 +1,188 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +# Copyright 2025 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. + +# Viscous-inviscid data structures +# Paul Dechamps + +import numpy as np +class Data: + """Data structure for quantities transfer + + Attributes: + ---------- + _name : str + Name of the data structure. + _type : str + Type of the data structure. + _ndim : int + Number of dimensions. + _nNodes : int + Number of nodes. + _nElms : int + Number of elements. + _nodesCoord : np.ndarray(float) + Coordinates of the nodes. + _elemsCoord : np.ndarray(float) + Coordinates of the elements. + _V : np.ndarray(float) + Velocity at the nodes. + _M : np.ndarray(float) + Mach number at the nodes. + _rho : np.ndarray(float) + Density at the nodes. + _connectNodes : np.ndarray(int) + Connect list of the nodes. + _connectElems : np.ndarray(int) + Connect list of the elements. + _nodeRows : np.ndarray(int) + Nodes rows. + _blowingVel : np.ndarray(float) + Blowing velocity at the elements. + init : bool + Initialization flag. + """ + def __init__(self, name:str, type:str, ndim:int): + self._name = name + self._type = type + self._ndim = ndim + + self._nNodes = 0 + self._nElms = 0 + + self.init = False + + # Setters + def updateVariables(self, M:np.ndarray, V:np.ndarray, Rho:np.ndarray): + """Update the variables of the data structure. + """ + if M.shape != self._M.shape: + raise Exception(f"Expected M with shape {self._M.shape}, got {M.shape}") + if V.shape != self._V.shape: + raise Exception(f"Expected V with shape {self._V.shape}, got {V.shape}") + if Rho.shape != self._rho.shape: + raise Exception(f"Expected Rho with shape {self._rho.shape}, got {Rho.shape}") + self._M[:] = M + self._rho[:] = Rho + self._V[:] = V + + def updateMesh(self, nodes:np.ndarray, elems:np.ndarray): + """Same as setMesh but performs several checks. + """ + if nodes.shape[0] != self._nNodes: + raise ValueError(f"updateMesh: Inconsistent number of nodes {nodes.shape[0]} vs {self._nNodes}") + if elems.shape[0] != self._nElms: + raise ValueError(f"updateMesh: Inconsistent number of elements {elems.shape[0]} vs {self._nElms}") + self._nodesCoord = nodes + self._elemsCoord = elems + self._broadcastInit() + + def setBlowingVelocity(self, blowingVel:np.ndarray): + if blowingVel.shape[0] != self._blowingVel.shape[0]: + raise Exception(f"Expected blowingVel with shape {self._blowingVel.shape[0]}, got {blowingVel.shape[0]}") + if blowingVel.shape[0] != self._nElms: + raise Exception(f"Expected blowingVel with shape {self._nElms}, got {blowingVel.shape}") + self._blowingVel[:] = blowingVel + + def setConnectLists(self, connectNodes:np.ndarray, connectElems:np.ndarray): + """Set the connect lists of the inviscid data. + """ + if connectNodes.shape[0] != self._nodesCoord.shape[0]: + raise Exception(f"Expected connectNodes with shape {self._nodesCoord.shape[0]}, got {connectNodes.shape[0]}") + if connectElems.shape[0] != self._elemsCoord.shape[0]: + raise Exception(f"Expected connectElems with shape {self._elemsCoord.shape[0]}, got {connectElems.shape[0]}") + self._connectNodes = np.array(connectNodes, dtype=int, order='C') + self._connectElems = np.array(connectElems, dtype=int, order='C') + + def setRows(self, nodesRow:np.ndarray, elemsRow:np.ndarray): + """Set the nodes rows of the inviscid data. + """ + if nodesRow.shape[0] != self._nNodes: + raise Exception(f"Expected nodesRow with shape {self._nNodes}, got {nodesRow.shape[0]}") + if elemsRow.shape[0] != self._nElms: + raise Exception(f"Expected elemsRow with shape {self._nElms}, got {elemsRow.shape[0]}") + self._nodeRows = np.array(nodesRow, dtype=int, order='C') + self._elemRows = np.array(elemsRow, dtype=int, order='C') + + def _setNodes(self, nodes:np.ndarray): + if nodes.shape[1] != 3: + raise Exception(f"Expected nodes with 3 dimensions, got {nodes.shape[1]}") + self._nodesCoord = np.array(nodes, dtype=float, order='C') + self._nNodes = nodes.shape[0] + + self._V = np.empty((self._nNodes, self._ndim), dtype=float, order='C') + self._M = np.empty(self._nNodes, dtype=float, order='C') + self._rho = np.empty(self._nNodes, dtype=float, order='C') + + def _setElems(self, elems:np.ndarray): + if elems.shape[1] != 3: + raise Exception(f"Expected elements with 3 dimensions, got {elems.shape[1]}") + self._elemsCoord = np.array(elems, dtype=float, order='C') + self._nElms = elems.shape[0] + + # Initialize to zero for the first iteration. Do not put empty. + self._blowingVel = np.zeros(self._nElms, dtype=float, order='C') + + # Getters + def getConnectListNodes(self): + return self._connectNodes + + def getBlowingVelocity(self): + return self._blowingVel + + def isinit(self): + return self.init + + def getType(self): + return self._type + + def getName(self): + return self._name + + def getnDim(self): + return self._ndim + + def getnNodes(self): + return self._nNodes + + def getnElms(self): + return self._nElms + + # Private functions + def _broadcastInit(self): + self.init = True + + # Functions to be implemented in the child classes + def getNodes(self): + raise NotImplementedError(f"getNodes not implemented for {self._name} of type {self._type}") + + def getElms(self): + raise NotImplementedError(f"getElms not implemented for {self._name} of type {self._type}") + + def getMach(self, reg:str): + raise NotImplementedError(f"getMach not implemented for {self._name} of type {self._type}") + + def getVelocity(self, reg:str): + raise NotImplementedError(f"getVelocity not implemented for {self._name} of type {self._type}") + + def getDensity(self, reg:str): + raise NotImplementedError(f"getDensity not implemented for {self._name} of type {self._type}") + + def computeStagPoint(self): + raise NotImplementedError(f"computeStagPoint not implemented for {self._name} of type {self._type}") + + def connectBlowingVel(self): + raise NotImplementedError(f"connectBlowingVel not implemented for {self._name} of type {self._type}") diff --git a/blast/interfaces/blDataInviscid.py b/blast/interfaces/blDataInviscid.py new file mode 100644 index 0000000000000000000000000000000000000000..0352031dfa8e63b8faf7df29679d7b5125eb91e1 --- /dev/null +++ b/blast/interfaces/blDataInviscid.py @@ -0,0 +1,233 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +# Copyright 2025 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. + +# Inviscid data structures +# Paul Dechamps + +import numpy as np +from blast.interfaces.blData import Data +class InviscidData(Data): + """ + Inviscid data structures + + Attributes + ---------- + _upNodes : np.ndarray + Array of the upper nodes index (in nodesCoord (for already sorted arrays if applicable)). + _lwNodes : np.ndarray + Array of the lower nodes index (in nodesCoord (for already sorted arrays if applicable)). + _upElems : np.ndarray + Array of the upper elements index (in elemsCoord (for already sorted arrays if applicable)). + _lwElems : np.ndarray + Array of the lower elements index (in elemsCoord (for already sorted arrays if applicable)). + """ + def __init__(self, _name, _type, _nDim): + Data.__init__(self, _name, _type, _nDim) + self.sidesIdentified = False + + def setMesh(self, nodes:np.ndarray, elems:np.ndarray): + """Set the mesh of the inviscid data. + + Parameters + ---------- + nodes : np.ndarray + The nodes coordinates. + elems : np.ndarray + The elements coordinates. + + Notes + ----- + The main difference with viscous data is that here we cannot + compute the elements coordinates automatically. + """ + self._setNodes(nodes) + self._setElems(elems) + self._broadcastInit() + + def getNodeRows(self): + """Get the nodes rows of the inviscid data. + """ + return self._nodeRows + + def getElmRows(self): + """Get the elements rows of the inviscid data. + """ + return self._elemRows + + def setUpperLower(self, upNodes:list, lwNodes:list, upElems:list, lwElems:list): + self._upNodes = np.array(upNodes, dtype=int, order='C') + self._lwNodes = np.array(lwNodes, dtype=int, order='C') + self._upElems = np.array(upElems, dtype=int, order='C') + self._lwElems = np.array(lwElems, dtype=int, order='C') + self.sidesIdentified = True + + def getNodes(self, reg:str='all'): + """Get the nodes coordinates. + + Parameters + ---------- + reg : str + The region of the nodes to get. + + Returns + ------- + np.ndarray + The nodes coordinates. + """ + if reg in ['all', 'fuselage', 'wake']: + return self._nodesCoord + elif reg == 'upper_fromLE': + #return self._nodesCoord[np.isin(self._nodeRows, self._upNodes), :] + nodes = np.empty((self._upNodes.shape[0], self._ndim), dtype=float, order='C') + for i in range(self._upNodes.shape[0]): + nodes[i, :] = self._nodesCoord[self._nodeRows == self._upNodes[i],:self._ndim] + return nodes + elif reg == 'lower_fromLE': + #return self._nodesCoord[np.isin(self._nodeRows, self._lwNodes), :] + nodes = np.empty((self._lwNodes.shape[0], self._ndim), dtype=float, order='C') + for i in range(self._lwNodes.shape[0]): + nodes[i, :] = self._nodesCoord[self._nodeRows == self._lwNodes[i],:self._ndim] + return nodes + else: + raise ValueError(f"getNodes: Unknown region {reg} for {self._name} of type {self._type}") + + def getElms(self, reg:str='all'): + """Get the elements coordinates. + + Parameters + ---------- + reg : str + The region of the elements to get. + + Returns + ------- + np.ndarray + The elements coordinates. + """ + if reg in ['all', 'fuselage', 'wake']: + return self._elemsCoord + elif reg == 'upper_fromLE': + #return self._elemsCoord[np.isin(self._elemRows, self._upElems), :] + elems = np.empty((self._upElems.shape[0], self._ndim), dtype=float, order='C') + for i in range(self._upElems.shape[0]): + elems[i, :] = self._elemsCoord[self._elemRows == self._upElems[i],:self._ndim] + return elems + elif reg == 'lower_fromLE': + #return self._elemsCoord[np.isin(self._elemRows, self._lwElems), :] + elems = np.empty((self._lwElems.shape[0], self._ndim), dtype=float, order='C') + for i in range(self._lwElems.shape[0]): + elems[i, :] = self._elemsCoord[self._elemRows == self._lwElems[i],:self._ndim] + return elems + else: + raise ValueError(f"getElms: Unknown region {reg} for {self._name} of type {self._type}") + + def getVelocity(self, reg:str='all'): + """Get the velocity. + + Parameters + ---------- + reg : str + The region of the velocity to get. + + Returns + ------- + np.ndarray + The velocity. + """ + if reg in ['all', 'fuselage', 'wake']: + return self._V + elif reg == 'upper_fromLE': + #return self._V[np.isin(self._nodeRows, self._upNodes), :] + v = np.empty((self._upNodes.shape[0], self._ndim), dtype=float, order='C') + for i, row in enumerate(self._upNodes): + for idim in range(self._ndim): + v[i, idim] = self._V[self._nodeRows == row, idim][0] + return v + elif reg == 'lower_fromLE': + #return self._V[np.isin(self._nodeRows, self._lwNodes), :] + v = np.empty((self._lwNodes.shape[0], self._ndim), dtype=float, order='C') + for i, row in enumerate(self._lwNodes): + for idim in range(self._ndim): + v[i, idim] = self._V[self._nodeRows == row, idim][0] + return v + else: + raise ValueError(f"getVelocity: Unknown region {reg} for {self._name} of type {self._type}") + + def getMach(self, reg:str='all'): + """Get the Mach number. + + Parameters + ---------- + reg : str + The region of the Mach number to get. + + Returns + ------- + np.ndarray + The Mach number. + """ + if reg in ['all', 'fuselage', 'wake']: + return self._M + elif reg == 'upper_fromLE': + #return self._M[np.isin(self._nodeRows, self._upNodes)] + m = np.empty(self._upNodes.shape[0], dtype=float) + for i, row in enumerate(self._upNodes): + m[i] = self._M[self._nodeRows == row][0] + return m + elif reg == 'lower_fromLE': + #return self._M[np.isin(self._nodeRows, self._lwNodes)] + m = np.empty(self._lwNodes.shape[0], dtype=float) + for i, row in enumerate(self._lwNodes): + m[i] = self._M[self._nodeRows == row][0] + return m + else: + raise ValueError(f"getMach: Unknown region {reg} for {self._name} of type {self._type}") + + def getDensity(self, reg:str='all'): + """Get the density. + + Parameters + ---------- + reg : str + The region of the density to get. + + Returns + ------- + np.ndarray + The density. + """ + if reg in ['all', 'fuselage', 'wake']: + return self._rho + elif reg == 'upper_fromLE': + #return self._rho[np.isin(self._nodeRows, self._upNodes)] + rho = np.empty(self._upNodes.shape[0], dtype=float) + for i, row in enumerate(self._upNodes): + rho[i] = self._rho[self._nodeRows == row][0] + return rho + elif reg == 'lower_fromLE': + #return self._rho[np.isin(self._nodeRows, self._lwNodes)] + rho = np.empty(self._lwNodes.shape[0], dtype=float) + for i, row in enumerate(self._lwNodes): + rho[i] = self._rho[self._nodeRows == row][0] + return rho + else: + raise ValueError(f"getDensity: Unknown region {reg} for {self._name} of type {self._type}") + + def connectBlowingVel(self): + """Connect the blowing velocity. + """ + self._blowingVel = self._blowingVel[self._connectElems.argsort()] diff --git a/blast/interfaces/blDataStructure.py b/blast/interfaces/blDataStructure.py deleted file mode 100644 index 7d9882bd1486a884e7baac3dd1b751b4c4cc1640..0000000000000000000000000000000000000000 --- a/blast/interfaces/blDataStructure.py +++ /dev/null @@ -1,398 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -# Copyright 2025 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. - -# Viscous-inviscid data structures -# Paul Dechamps - -import numpy as np -class Group: - def __init__(self, _name, _type, _nDim): - """Initialize the group. - - args: - ---- - _name: str - Name of the group - _nDim: int - Number of dimensions - """ - self.name = _name - self.type = _type - self.nDim = _nDim - self.stagPoint = None - self.sidesIdentified = False - - self.upNodes = None - self.lwNodes = None - self.upElems = None - self.lwElems = None - - def initStructures(self, nPoints, nElems): - """Initialize the data structures. - - args: - ---- - nPoints: int - Number of points in the region - """ - self.nPoints = nPoints - self.nElems = nElems - self.nodesCoord = np.zeros((self.nPoints,4)) - self.elemsCoord = np.zeros((self.nElems,4)) - - self.V = np.zeros((self.nPoints,3)) - self.M = np.zeros(self.nPoints) - self.Rho = np.zeros(self.nPoints) - self.blowingVel = np.zeros(self.nElems) - - def setUpperLower(self, _upNodes, _lowNodes, _upElems, _lwElems): - """Set the upper and lower nodes and elements for the group. - - args: - ---- - _upNodes: np.ndarray - Array of upper node indices - _lowNodes: np.ndarray - Array of lower node indices - _upElems: np.ndarray - Array of upper element indices - _lwElems: np.ndarray - Array of lower element indices - """ - self.upNodes = _upNodes - self.lwNodes = _lowNodes - self.upElems = _upElems - self.lwElems = _lwElems - - def setStagPoint(self, _stagPoint): - """Set the stagnation point. - - args: - ---- - _stagPoint: int - Index of the stagnation point - """ - self.stagPoint = _stagPoint - - ##### NODES GETTERS ##### - def getUpperNodes(self): - """Return the upper nodes. - """ - if self.name == 'iWing': - nodes = np.empty((self.upNodes.shape[0], self.nDim), dtype=float, order='C') - for i in range(self.upNodes.shape[0]): - nodes[i, :] = self.nodesCoord[self.nodesCoord[:,-1] == self.upNodes[i],:self.nDim] - return nodes - elif 'vAirfoil' in self.name: - return np.array(self.nodesCoord[:np.argmin(self.nodesCoord[:,0])+1, :self.nDim], dtype=float) - else: - raise RuntimeError(f'Lower nodes can only be retrieved for iWing of vAirfoil, not {self.name}') - - def getLowerNodes(self): - """Return the lower nodes. - """ - if self.name == 'iWing': - nodes = np.empty((self.lwNodes.shape[0], self.nDim), dtype=float, order='C') - for i in range(self.lwNodes.shape[0]): - nodes[i, :] = self.nodesCoord[self.nodesCoord[:,-1] == self.lwNodes[i],:self.nDim] - return nodes - elif 'vAirfoil' in self.name: - return np.array(self.nodesCoord[np.argmin(self.nodesCoord[:,0])+1:, :self.nDim], dtype=float) - else: - raise RuntimeError(f'Lower nodes can only be retrieved for iWing of vAirfoil, not {self.name}') - - ##### ELEMENTS GETTERS ##### - def getUpperElems(self): - """Return the upper elements. - """ - if self.name == 'iWing': - elems = np.empty((self.upElems.shape[0], self.nDim), dtype=float, order='C') - for i in range(self.upElems.shape[0]): - elems[i, :] = self.elemsCoord[self.elemsCoord[:,-1] == self.upElems[i],:self.nDim] - return elems - elif 'vAirfoil' in self.name: - return np.array(self.elemsCoord[:np.argmin(self.nodesCoord[:,0]), :self.nDim], dtype=float) - else: - raise RuntimeError(f'Upper elements can only be retrieved for iWing or vAirfoil, not {self.name}') - - - def getLowerElems(self): - """Return the lower elements. - """ - if self.name == 'iWing': - elems = np.empty((self.lwElems.shape[0], self.nDim), dtype=float, order='C') - for i in range(self.lwElems.shape[0]): - elems[i, :] = self.elemsCoord[self.elemsCoord[:,-1] == self.lwElems[i],:self.nDim] - return elems - elif 'vAirfoil' in self.name: - return np.array(self.elemsCoord[np.argmin(self.nodesCoord[:,0]):, :self.nDim], dtype=float) - else: - raise RuntimeError(f'Lower elements can only be retrieved for iWing of vAirfoil, not {self.name}') - - ##### VELOCITY GETTERS ##### - def getUpperVelocity(self): - if self.name == 'iWing': - v = np.empty((self.upNodes.shape[0], self.nDim), dtype=float, order='C') - for i, row in enumerate(self.upNodes): - for idim in range(self.nDim): - v[i, idim] = self.V[self.nodesCoord[:,3] == row, idim][0] - return v - else: - raise RuntimeError(f'Upper velocity can only be retrieved for iWing, not {self.name}') - - def getLowerVelocity(self): - if self.name == 'iWing': - v = np.empty((self.lwNodes.shape[0], self.nDim), dtype=float, order='C') - for i, row in enumerate(self.lwNodes): - for idim in range(self.nDim): - v[i, idim] = self.V[self.nodesCoord[:,3] == row, idim][0] - return v - else: - raise RuntimeError(f'Lower velocity can only be retrieved for iWing, not {self.name}') - - ##### MACH NUMBER GETTERS ##### - def getUpperMach(self): - if self.name == 'iWing': - m = np.empty(self.upNodes.shape[0], dtype=float) - for i, row in enumerate(self.upNodes): - m[i] = self.M[self.nodesCoord[:,3] == row][0] - return m - else: - raise RuntimeError(f'Upper Mach number can only be retrieved for iWing, not {self.name}') - - def getLowerMach(self): - if self.name == 'iWing': - m = np.empty(self.lwNodes.shape[0], dtype=float) - for i, row in enumerate(self.lwNodes): - m[i] = self.M[self.nodesCoord[:,3] == row][0] - return m - else: - raise RuntimeError(f'Lower Mach number can only be retrieved for iWing, not {self.name}') - - ##### DENSITY GETTERS ##### - def getUpperDensity(self): - if self.name == 'iWing': - rho = np.empty(self.upNodes.shape[0], dtype=float) - for i, row in enumerate(self.upNodes): - rho[i] = self.Rho[self.nodesCoord[:,3] == row][0] - return rho - else: - raise RuntimeError(f'Upper density can only be retrieved for iWing, not {self.name}') - - def getLowerDensity(self): - if self.name == 'iWing': - rho = np.empty(self.lwNodes.shape[0], dtype=float) - for i, row in enumerate(self.lwNodes): - rho[i] = self.Rho[self.nodesCoord[:,3] == row][0] - return rho - else: - raise RuntimeError(f'Lower density can only be retrieved for iWing, not {self.name}') - - ##### BLOWING VELOCITY GETTERS ##### - def getUpperBlowingVel(self): - if self.name == 'iWing': - v = np.empty(self.upElems.shape[0], dtype=float) - for i, row in enumerate(self.upElems): - v[i] = self.blowingVel[self.elemsCoord[:,3] == row][0] - return v - elif 'vAirfoil' in self.name: - return np.array(self.blowingVel[:np.argmin(self.nodesCoord[:,0])], dtype=float) - else: - raise RuntimeError(f'Upper blowing velocity can only be retrieved for iWing, not {self.name}') - - def getLowerBlowingVel(self): - if self.name == 'iWing': - v = np.empty(self.lwElems.shape[0], dtype=float) - for i, row in enumerate(self.lwElems): - v[i] = self.blowingVel[self.elemsCoord[:,3] == row][0] - return v - elif 'vAirfoil' in self.name: - return np.array(self.blowingVel[np.argmin(self.nodesCoord[:,0]):], dtype=float) - else: - raise RuntimeError(f'Lower blowing velocity can only be retrieved for iWing, not {self.name}') - - def updateVars(self, _V, _M, _Rho): - """Update the velocity, Mach number and density of the boundary layer. - - args: - ---- - _V: np.ndarray - Velocity - _M: np.ndarray - Mach number - _Rho: np.ndarray - Density - """ - self.V = _V - self.M = _M - self.Rho = _Rho - - def getNodesCoord(self, reg:str): - """Return the coordinates of the nodes in the direction dim - - args: - ---- - reg: str - Region of the boundary layer (upper, lower, wake) - dim: int - Direction of the coordinates (x, y, z) - """ - if reg in ["wake", "fuselage"]: - return self.nodesCoord - elif reg == "upper": - return np.flip(self.nodesCoord[:self.stagPoint+1,:], axis=0) - elif reg == "lower": - return self.nodesCoord[self.stagPoint:,:] - else: - raise RuntimeError('Invalid region', reg) - - def getNodesRow(self, reg:str): - """Return the row of the nodes - - args: - ---- - reg: str - Region of the boundary layer (upper, lower, wake) - """ - if reg in ["wake", "fuselage"]: - return self.connectListNodes - elif reg == "upper": - return np.flip(self.connectListNodes[:self.stagPoint+1]) - elif reg == "lower": - return self.connectListNodes[self.stagPoint:] - else: - raise RuntimeError('Invalid region', reg) - - def getVelocity(self, reg:str, dim:int): - """Return the velocity in the direction dim - - args: - ---- - reg: str - Region of the boundary layer (upper, lower, wake) - dim: int - Direction of the velocity - """ - if reg in ["wake", "fuselage"]: - return self.V[:,dim] - elif reg == "upper": - return np.flip(self.V[:self.stagPoint+1,dim]) - elif reg == "lower": - return self.V[self.stagPoint:,dim] - else: - raise RuntimeError('Invalid region', reg) - - def getVt(self, reg:str): - """Return the velocity in the direction tangential - to the boundary layer - - args: - ---- - reg: str - Region of the boundary layer (upper, lower, wake) - """ - vt = np.sqrt(self.V[:,0]**2 + self.V[:,self.nDim-1]**2) - if reg in ["wake", "fuselage"]: - return vt - elif reg == "upper": - return np.flip(vt[:self.stagPoint+1]) - elif reg == "lower": - return vt[self.stagPoint:] - else: - raise RuntimeError('Invalid region', reg) - - def getMach(self, reg:str): - """Return the Mach number at the edge of the boundary layer - - args: - ---- - reg: str - Region of the boundary layer (upper, lower, wake) - """ - if reg in ["wake", "fuselage"]: - return self.M - elif reg == "upper": - return np.flip(self.M[:self.stagPoint+1]) - elif reg == "lower": - return self.M[self.stagPoint:] - else: - raise RuntimeError('Invalid region', reg) - - def getRho(self, reg:str): - """Return the density at the edge of the boundary layer - - args: - ---- - reg: str - Region of the boundary layer (upper, lower, wake) - """ - if reg in ["wake", "fuselage"]: - return self.Rho - elif reg == "upper": - return np.flip(self.Rho[:self.stagPoint+1]) - elif reg == "lower": - return self.Rho[self.stagPoint:] - else: - raise RuntimeError('Invalid region', reg) - - def getConnectList(self, reg:str): - """Return the connectivity list of the elements - - args: - ---- - reg: str - Region of the boundary layer (upper, lower, wake) - """ - if reg in ["wake", "fuselage"]: - return self.connectListElems - elif reg == "upper": - return np.flip(self.connectListElems[:self.stagPoint]) - elif reg == "lower": - return self.connectListElems[self.stagPoint:] - else: - raise RuntimeError('Invalid region', reg) - - def setConnectList(self, _connectListNodes, _connectListElems): - """Set connectivity lists for viscous structures. - - args: - ---- - _connectListNodes: np.ndarray - Connectivity list for nodes - _connectListElems: np.ndarray - Connectivity list for elements - """ - if self.name != 'iWing' and self.name != 'iWake': - raise RuntimeError('Cannot set connectivity lists for viscous structures', self.name) - self.connectListNodes = _connectListNodes - self.connectListElems = _connectListElems - - def connectBlowingVel(self): - """Connect blowing velocities for viscous structures - using the connectivity list. - """ - if self.name != 'iWing' and self.name != 'iWake': - raise RuntimeWarning('Can not connect blowing velocities for viscous structure. Maybe it was not set.', self.name) - self.blowingVel = self.blowingVel[self.connectListElems.argsort()] - - def computeStagPoint(self): - """Compute the stagnation point of the boundary layer. - """ - if 'iWake' in self.name or 'vWake' in self.name or self.type == 'fuselage': - self.stagPoint = None - else: - self.stagPoint = np.argmin(np.sqrt(self.V[:,0]**2 + self.V[:,self.nDim-1]**2)) diff --git a/blast/interfaces/blDataViscous.py b/blast/interfaces/blDataViscous.py new file mode 100644 index 0000000000000000000000000000000000000000..53d5ee35ce2056446146d385481697e9e05d610e --- /dev/null +++ b/blast/interfaces/blDataViscous.py @@ -0,0 +1,229 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +# Copyright 2025 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. + +# Viscous data structures +# Paul Dechamps + +import numpy as np +from blast.interfaces.blData import Data +class ViscousData(Data): + def __init__(self, _name, _type, _nDim): + Data.__init__(self, _name, _type, _nDim) + self._stagPoint = None + self._leNodeIndex = None + + def setMesh(self, nodes:np.ndarray, elems:np.ndarray|None = None): + self._setNodes(nodes) + self._leNodeIndex = np.argmin(self._nodesCoord[:, 0]) + + if elems is None: + n_nodes = self.getnNodes() + elems = np.empty((n_nodes - 1, 3), dtype=float, order='C') + elems[:, :] = 0.5 * (self._nodesCoord[:-1, :] + self._nodesCoord[1:, :]) + self._setElems(elems) + + def getNodes(self, reg:str='all')->np.ndarray: + """Return the nodes coordinates + + args: + ---- + reg: str + Region of the boundary layer (all, upper, lower, wake), default is all. + """ + if reg in ["all", "wake", "fuselage"]: + return self._nodesCoord + elif reg == "upper": + return np.flip(self._nodesCoord[:self._stagPoint+1,:], axis=0) + elif reg == "lower": + return self._nodesCoord[self._stagPoint:,:] + elif reg == "upper_fromLE": + return self._nodesCoord[:self._leNodeIndex+1,:] + elif reg == "lower_fromLE": + return self._nodesCoord[self._leNodeIndex+1:,:] + else: + raise RuntimeError(f'Unrecognized region {reg} for {self.getName()} of type {self.getType()}') + + def getElms(self, reg:str='all')->np.ndarray: + """Return the elements coordinates + + args: + ---- + reg: str + Region of the boundary layer (all, upper, lower, wake), default is all. + """ + if reg in ["all", "wake", "fuselage"]: + return self._elemsCoord + elif reg == "upper": + return np.flip(self._elemsCoord[:self._stagPoint,:], axis=0) + elif reg == "lower": + return self._elemsCoord[self._stagPoint:,:] + elif reg == "upper_fromLE": + return self._elemsCoord[:self._leNodeIndex,:] + elif reg == "lower_fromLE": + return self._elemsCoord[self._leNodeIndex:,:] + else: + raise RuntimeError(f'Unrecognized region {reg} for {self.getName()} of type {self.getType()}') + + def getVelocity(self, reg:str='all', dim:int | list = [0,1,2]): + """Return the velocity in the direction dim + + args: + ---- + reg: str + Region of the boundary layer (all, upper, lower, wake), default is all. + dim: int + Direction of the velocity, default is [0,1,2] + """ + if reg in ["all", "wake", "fuselage"]: + return self._V[:,dim] + elif reg == "upper": + return np.flip(self._V[:self._stagPoint+1,dim]) + elif reg == "lower": + return self._V[self._stagPoint:,dim] + # elif reg == "upper_fromLE": + # return np.flip(self._V[:self._leNodeIndex,dim]) + # elif reg == "lower_fromLE": + # return self._V[self._leNodeIndex:,dim] + else: + raise RuntimeError(f'Unrecognized region {reg} for {self.getName()} of type {self.getType()}') + + # Boundary layer quantities + def getVt(self, reg:str="all"): + """Return the velocity in the direction tangential + to the boundary layer + + args: + ---- + reg: str + Region of the boundary layer (upper, lower, wake) + """ + vt = np.sqrt(self._V[:,0]**2 + self._V[:,self._ndim-1]**2) + if reg in ["all", "wake", "fuselage"]: + return vt + elif reg == "upper": + return np.flip(vt[:self._stagPoint+1]) + elif reg == "lower": + return vt[self._stagPoint:] + else: + raise RuntimeError(f'Unrecognized region {reg} for {self.getName()} of type {self.getType()}') + + def getMach(self, reg:str="all"): + """Return the Mach number + """ + if reg in ["all", "wake", "fuselage"]: + return self._M + elif reg == "upper": + return np.flip(self._M[:self._stagPoint+1]) + elif reg == "lower": + return self._M[self._stagPoint:] + else: + raise RuntimeError(f'Unrecognized region {reg} for {self.getName()} of type {self.getType()}') + + def getBlowingVelocity(self, reg:str="all"): + """Return the blowing velocity + """ + if reg in ["all", "wake", "fuselage"]: + return self._blowingVel + elif reg == "upper": + return np.flip(self._blowingVel[:self._stagPoint]) + elif reg == "lower": + return self._blowingVel[self._stagPoint:] + elif reg == "upper_fromLE": + return self._blowingVel[:self._leNodeIndex] + elif reg == "lower_fromLE": + return self._blowingVel[self._leNodeIndex:] + else: + raise RuntimeError(f'Unrecognized region {reg} for {self.getName()} of type {self.getType()}') + + def getDensity(self, reg:str="all"): + """Return the density + """ + if reg in ["all", "wake", "fuselage"]: + return self._rho + elif reg == "upper": + return np.flip(self._rho[:self._stagPoint+1]) + elif reg == "lower": + return self._rho[self._stagPoint:] + else: + raise RuntimeError(f'Unrecognized region {reg} for {self.getName()} of type {self.getType()}') + + def getNodeRows(self, reg:str='all'): + """Get the nodes rows of the inviscid data. + """ + if reg in ["all", "wake", "fuselage"]: + return self._nodeRows + elif reg == 'upper': + return np.flip(self._nodeRows[:self._stagPoint+1], axis=0) + elif reg == 'lower': + return self._nodeRows[self._stagPoint:] + else: + raise ValueError(f"Unknown region {reg} for {self.getName()} of type {self.getType()}") + + def getElmsRows(self, reg:str='all'): + """Get the elements rows of the inviscid data. + """ + if reg in ["all", "wake", "fuselage"]: + return self._elmsRows + elif reg == 'upper': + return np.flip(self._elmsRows[:self._stagPoint], axis=0) + elif reg == 'lower': + return self._elmsRows[self._stagPoint:] + else: + raise ValueError(f"Unknown region {reg} for {self.getName()} of type {self.getType()}") + + def getNodeConnectList(self, reg:str='all'): + """Get the connect lists of the inviscid data. + """ + if reg in ["all", "wake", "fuselage"]: + return self._connectNodes + elif reg == 'upper': + return np.flip(self._connectNodes[:self._stagPoint+1], axis=0) + elif reg == 'lower': + return self._connectNodes[self._stagPoint:] + else: + raise ValueError(f"Unknown region {reg} for {self.getName()} of type {self.getType()}") + + def getElmsConnectList(self, reg:str='all'): + """Get the connect lists of the inviscid data. + """ + if reg in ["all", "wake", "fuselage"]: + return self._connectElems + elif reg == 'upper': + return np.flip(self._connectElems[:self._stagPoint], axis=0) + elif reg == 'lower': + return self._connectElems[self._stagPoint:] + else: + raise ValueError(f"Unknown region {reg} for {self.getName()} of type {self.getType()}") + + # Stagnation point + def getStagPoint(self): + return self._stagPoint + + def setStagPoint(self, _stagPoint): + if _stagPoint is not None and _stagPoint > self.getnNodes(): + raise ValueError("Stagnation point index out of range") + self._stagPoint = _stagPoint + + def computeStagPoint(self): + """Compute the stagnation point of the boundary layer. + """ + if 'iWake' in self.getName()\ + or 'vWake' in self.getName()\ + or self.getType() == 'fuselage': + self.setStagPoint(None) + else: + self.setStagPoint(np.argmin(self.getVt('all')))